TRANSONIC  THERMAL  BLOOMING 


Edwin  Fenton  Carey 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOG 
MONTEREY,  CALIFORNIA  9394(1 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


THESIS 

TRANSONIC  THERMAL  BLOOMING 

by 

Edwin  Fenton  Carey,   Jr. 

March  1976 

Thesis 

Advisor:                     A.  E. 

Fuhs 

Approved  for  public  release;  distribution  unlimited. 


T173123 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF   THIS  PACE  (When  Dete  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1.     REPORT  NUMBER 

2.  GOVT  ACCESSION  NO. 

1.     RECIPIENT'S  CATALOG  NUMBER 

4.     TITLE  (and  Subttttm) 

Transonic  Thermal  Blooming 

5.     TYPE  OF   REPORT  A  PERIOD  COVERED 

Ph.D.    Thesis 
March   1976 

•  ■    PERFORMING  ORG.   REPORT  NUMBER 

7.     AUTHORS 

Edwin  Fenton  Carey,    Jr. 

■  .     CONTRACT  OR  GRANT  NUMBERS 

9.     PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Naval   Postgraduate   School 
Monterey,    California      9  3940 

!0.     PROGRAM  ELEMENT.  PROJECT,   TASK 
AREA  *  WORK  UNIT  NUMBERS 

I  1.     CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval   Postgraduate   School 
Monterey,    California      93940 

12.     REPORT   DATE 

March   1976 

13.     NUMBER  OF  PAGES 
122 

14.     MONITORING  AGENCY  NAME  *    ADORES^//  dlllarent  from  Controlling  Oil  Ice') 

IS.    SECURITY  CLASS,  (ol  thta  riporl) 

Unclassified 

1Sa.     OECLASSIFI CATION/ DOWNGRADING 
SCHEDULE 

1«.     DISTRIBUTION   STATEMENT  (of  (hit  Riper!) 

Approved   for  public   release;    distribution   unlimited. 

17.     DISTRIBUTION  STATEMENT  (el  the  abatrect  entered  In  Block  30,  il  dltfmrent  from  Report) 

18.     SUPPLEMENTARY  NOTES 

19.     KEY  WORDS  (Continue  on  rererae  aide  II  nacaatery  and  Identity  by  block  number) 

Transonic   Thermal  Blooming                             High  Energy  Laser 
Subsonic   Thermal   Blooming 
Supersonic  Thermal   Blooming 
Atmospheric  Propagation  of   Laser   Beam 
Sonic   Flow  with   Heat  Addition 

20.     ABSTRACT  (Continue  on  rererae  aide  It  neceaaawy  end  Identity  by  block  number) 

According   to   the   linearized   solutions    for   thermal  blooming, 
the   density  perturbations   become   infinite    (i.e.    "catastrophic" 
defocusing)    as    the  Mach  number   approaches   unity.      However,    the 
nonlinearities    in   the   transonic   equations    cutoff   the   trend   to 
infinity,    and   the   values   of    the   flow  perturbation  quantities 
are   finite.      The   nonlinear   equations   with  heat  addition  are 
transformed   into   simple   linear   algebraic   equations    through   the 

dd  ,; 

(Page    1) 


FaT7,  1473 


EDITION  OF  1  NOV  •■  IS  OBSOLETE 

S/N    0  102-014- 6601    | 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAOE  (When  Dmta  Kntered) 


UNCLASSIFIED 


fLiCUWlTV    CL  ASSlFICATlQN    OF    THIS    P>GEf*>»n    Orlm   EnHnJ 


(20.   CONTINUED) 

specification  of  the  streamline  geometry  in  the  heat  release 

region.   At  a  Mach  number  of  unity,  streamtube  area  variation 

was  found  to  be  directly  proportional  to  the  change  in  total 

temperature.   A  steady,  two  -dimensional  mixed  flow  solution 

has  been  found  for  the  transonic  thermal  blooming  problem. 

The  solution  for  the  density  perturbations  within  a  laser 

beam  at  a  Mach  number  of  precisely  unity  is  given.   For  a 

7        2 
Gaussian  beam  with  an  intensity  of  3.333x10  Watts/m  and  an 

-7   -I 

atmospheric  absorption  of  8.0x10    cm   the  maximum  fractional 

density  perturbation  is  1.028x10   .   The  transonic  thermal 
blooming  problem  does  not  pose  as  serious  a  problem  as 
previously  anticipated. 


DD      Form       1473 

.  1  Jan  73  UNCLASSIFIED — 

S/N       0102-014-6601  SECURITY   CLASSIFICATION   OF   THIS  P»CEfW.n  Dmtm  Enfr.d) 


Transonic  Thermal  Blooming 

by 

Edwin  Fenton  Carey,  Jr. 

Lieutenant  Commander,  tJnited  States  Navy 

B.M.A.E.,  University  of  Delaware,  1967 

M.M.A.E.,  University  of  Delaware,  1970 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

DOCTOR  OF  PHILOSOPHY 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
March  1976 


Tk€i>5 


DUDLEY   KNOX  LIBRARY 
NAVAL  POSTGRADUATE  8CHOOQ 
MONTEREY,  CALIFORNIA  9394Q 


ABSTRACT 


According  to  the  linearized  solutions  for  thermal 
blooming,  the  density  perturbations  become  infinite  (i.e. 
"catastrophic"  defocusing)  as  the  Mach  number  approaches 
unity.   However,  the  nonlinear ities  in  the  transonic  equations 
cutoff  the  trend  to  infinity,  and  the  values  of  the  flow 
perturbation  quantities  are  finite.   The  nonlinear  equations 
with  heat  addition  are  transformed  into  simple  linear 
algebraic  equations  through  the  specification  of  the 
streamline  geometry  in  the  heat  release  region.   At  a  Mach 
number  of  unity,  streamtube  area  variation  was  found  to  be 
directly  proportional  to  the  change  in  total  temperature. 
A  steady,  two-dimensional  mixed  flow  solution  has  been  found 
for  the  transonic  thermal  blooming  problem.   The  solution 
for  the  density  perturbations  within  a  laser  beam  at  a  Mach 

number  of  precisely  unity  is  given.   For  a  Gaussian  beam 

7        2 
with  an  intensity  of  3.333x10   Watts/m  and  an  atmospheric 

-7    -1 
absorption  of  8.0x10    cm   the  maximum  fractional  density 

_  c 

perturbation  is  1.02  8x10   .   The  transonic  thermal  blooming 
problem  does  not  pose  as  serious  a  problem  as  previously 
anticipated. 
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NOTATION 

English 

A       area,  square  centimeters 

a       speed  of  sound,  meters/second 

b       function,  dimensionless 

C       line  contour  around  integration  region  R 

c       speed  of  light,  meters/second 

Cp      specific  heat  at  constant  pressure, 
Joules/kilogram  mass  -  degree  Kelvin 

Cv      specific  heat  at  constant  volume, 

Joules/kilogram  mass  -  degree  Kelvin 

f       function,  dimensionless 

G       Green's  function,  dimensionless;  also  specified 
function,  dimensionless 

H       Heat  quantity,  dimensionless 

h       specific  enthalpy,  Joules/kilogram  mass 

h(x,y)   heat  addition  function,  Watts/cubic  centimeters 

I       laser  beam  intensity,  Watts/square  centimeter; 
also  integral  function  representation 

I(x)     unit  function;  zero  for  x  <  0  and  unity  for  x   0 

k       Gladstone-Dale  constant,  cubic  meters/kilogram  mass 

L       scaling  factor,  dimensionless;  length,  meters; 

and  functional  representation  of  normalized  velocity, 
dimensionless 

I  characteristic  length  (subscript  1  for  x-direction 

and  2  for  y-direction) ,  centimeters 

M       Mach  number,  dimensionless 

N       Number  of  waves,  dimensionless 

n       normal  coordinate,  centimeters;  also  index  of 
refraction,  dimensionless 


p       static  pressure,  Newtons/square  meter 

Pr      Prandtl  number,  dimensionless 

Q       heat  source,  Watts/cubic  centimeter 

q       heat  source,  Watts/square  centimeter;  also, 

strength  of  line  heat  source,  Watts/centimeter 

R       universal  gas  constant,  Joules/kilogram  mass  - 

degrees  Kelvin  ;  radius  of  curvature,  centimeters; 
and  region  of  integration 

r       distance,  centimeters 

Re      Reynolds  number,  dimensionless 

S       line  contour  around  shock  waves;  also,  distance 
along  a  characteristic,  centimeters 

s       specific  entropy,  Joules/kilogram  mass  -  degree 
Kelvin  ;  also  streamline  coordinate 

T  static  temperature,  degrees  Kelvin 

t  maximum  streamtube  thickness,  centimeters 

U  velocity  vector,  meters/second 

U,u  speed  (x-direction) ,  meters/second 

v  speed  (y-direction) ,  meters/second 

w       laser  beam  spot  size,  centimeters; 
also  speed,  meters/second 

x       coordinate,  centimeters 

y       coordinate,  centimeters 

Z       streamtube  shape,  centimeters 


Greek 

a       atmospheric  absorption,  1/centimeters 

2  1/2 
3       variable,  dimensionless;  (1-M  )     for  subsonic 


flow,  dimensionless; 
and  (M2  -  1) 1/2   for  s 

ratio  of  specific  heats,  dimensionless 


2     1/2 
and  (M  -  1)  '       for  supersonic  flow,  dimensionless 
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A       small  change,  dimensionless 

<5        slope,  dimensionless;  also  delta  (impulse) function 

£       dimensionless  small  quantity,  0  <  e  <<  1 

n       absolute  viscosity,  kilogram  mass/meter-second; 
also  coordinate  (y-direction) ,  centimeters 

0       flow  angle,  radians 

k       thermal  conductivity,  Watts/centimeter  - 
degree  Kelvin 

X       wavelength,  centimeters 

y       Mach  angle  defined  by  arcsin(l/M) 

v       inverse  radius  of  curvature,  1/centimeter 

£       coordinate  (x-direction) ,  centimeters;  also 

transonic  similarity  parameter,  dimensionless 

p       density,  kilogram  mass/cubic  meter 

Z       summation 

a  standard  deviation  of  Gaussian  beam,  centimeters 

t       thickness  ratio,  dimensionless;  also,  vibrational 
relaxation  time,  seconds 

$  viscous  dissipation  function,  Watts/cubic  centimeter 

<j>  velocity  potential,  dimensionless 

i>  viscous  force  vector,  Newtons/cubic  meter 

ft  angular  velocity,  radians/second 

Superscripts 

(  ) '     prime,  dimensionless  quantity 

(  )      bar,  normalized  quantity;  also  average  value  of  a 
quantity 

(~)      tilde,  perturbation  quantity 
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Subscripts 

a,b     indicate  values  of  a  quantity  in  front  of  and 
behind  a  shock  wave 

i       x-direction  or  streamline  index  -  first  subscript  ■ 
(i  = 1  denotes  freestream  conditions) 

j       y-direction  or  normal  index  -  second  subscript  - 
( j  = 1  denotes  center line  conditions) 

L       linear  value 

0       total  or  stagnation  condition;  also,  peak  value  of 
an  intensity  distribution 

w       wall  or  surface  value 

»       freestream  condition 

1,2,     first,  second,  etc.  perturbation  quantities 


Miscellaneous 

V  Del  operator,  1/centimeters 

2 

V  Laplace  operator,  1/square  centimeters 
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I.   INTRODUCTION 

The  optical  quality  of  a  propagating  laser  beam  is  de- 
graded when  there  are  density  gradients  in  the  medium.   A 
laser  beam  propagating  through  an  absorbing  medium  creates 
local  changes  in  the  density  which  affect  the  index  of  refrac- 
tion of  the  medium.   Therefore,  the  density  gradients  will 
refract  the  light  into  regions  of  higher  index  of  refraction 
(higher  density)  causing  defocusing  of  the  beam.   This  pro- 
cess, which  is  known  as  thermal  blooming,  has  been  studied 
by  many  authors  [Refs.  1,  2,  and  3]  to  determine  its  effect 
on  a  laser  beam  propagating  through  various  media. 

When  a  laser  beam  is  slued,  a  complicated  coupled  pro- 
cess of  thermal  expansion  and  fluid  flow  (relative  winds) 
takes  place.   Depending  on  the  rate  of  sluing,  the  beam  will 
experience  subsonic,  transonic,  and  supersonic  winds  at  vari- 
ous radial  distances  from  the  laser  source.   The  subsonic 
and  supersonic  thermal  blooming  problems  have  been 
solved  [Refs.  4-8],   A  computer  program  (DIDER)  has  been 
developed  and  experimentally  verified  [Ref.  9]  for  beam 
interaction  with  supersonic  flows  internal  to  a  Gas  Dynamic 
Laser. 

A  solution  for  the  transonic  regime,  and,  in  particular, 
the  case  when  the  frees tream  Mach  number  is  precisely  unity, 
has  been  of  considerable  interest  because  the  addition  of 
heat  can  lead  to  extremely  strong  density  gradients.   As 
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with  the  theoretical  investigation  of  transonic  flows  over 
airfoils  and  wedges,  the  linear  small  perturbation  flow 
equations  valid  for  subsonic  flow  (elliptic  equations)  and 
supersonic  flow  (hyperbolic  equations)  are  no  longer  complete 
in  the  transonic  regime. 

Many  authors  [Refs.  10  -  15]  have  studied  various  methods 
of  solving  the  transonic  flow  equations  without  heat  addition. 
The  assumptions  based  on  known  transonic  behavior  about  air- 
foils and  wedges  are  used  to  simplify  the  governing  equations 
and,  further,  to  obtain  approximate  numerical  solutions. 
Often,  the  assumptions  greatly  restrict  the  types  of  transonic 
flows  that  can  be  solved.   For  example,  in  the  parabolic 
method  [Ref.  16],  the  transonic  flow  equation  is  transformed 
into  a  parabolic  differential  equation.   Since  the  solutions 
to  parabolic  differential  equations  are  well  established  in 
the  literature,  this  method  may  be  used  when  the  acceleration 
or  deceleration  is  approximately  constant  in  the  region  of 
interest  (i.e.  across  the  cord  of  many  airfoils) ;  however, 
this  method  is  not  valid  in  the  region  of  the  stagnation 
point  or  where  the  flow  acceleration  changes  sign. 

With  heat  addition,  the  basic  flow  equations  for  transonic 
flow  become  more  complicated.   Zierep  [Ref.  17]  has  solved 
the  case  when  the  Mach  number  is  precisely  unity;  but,  in 
solving  the  equation,  he  has  reduced  the  problem  to  that  of 
one  dimension.   Others  [Refs.  18  and  19]  have  solved  the 
nonsteady  transonic  equation  with  heat  addition;  however, 
the  results  are  valid  only  for  a  short  time  after  beam  turn 
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on  (i.e.  no  steady  state  solution) .   Even  for  the  short  time 
solutions,  extrapolation  techniques  between  the  property- 
values  at  the  critical  points  on  the  subsonic  side  and  the 
supersonic  side  are  employed  to  obtain  the  values  across 
the  transonic  regime.   Experimental  investigations  are  being 
conducted  [Ref.  20]  to  determine  the  influence  of  transonic 
sluing  on  thermal  blooming.   These  experiments  indicate  that 
a  steady  two-dimensional  solution  is  possible.   Further 
experimentation  will  provide  better  guidelines  for  a 
theoretical  analysis  of  the  problem. 

In  this  thesis,  transonic  thermal  blooming  in  steady 
two-dimensional  invisid  flow  is  formulated  and  solved  for  a 
Gaussian  heat  release  distribution.   In  particular,  it  examines 
the  problem  when  the  freestream  Mach  number  is  precisely 
unity.   This  is  accomplished  by  using  a  method  developed  by 
Broadbent  [Refs.  21  -  25]  in  which  the  streamlines  through- 
out the  flow  are  adjusted  in  such  a  manner  that  the  heat 
addition  necessary  to  achieve  these  streamlines  equals  the 
required  heat  distribution  and  boundary  conditions.   The 
method  is  exact  since  specifying  the  streamlines  reduces 
the  nonlinear  coupled  momentum  equations  in  natural  (stream- 
line) coordinates  to  simple  algebraic  linear  difference  equa- 
tions that  can  be  "marched"  through  the  entire  flow  field. 
Boundary  conditions  are  the  freestream  properties  upstream 
of  the  heat  release  region  and  those  on  a  bounding  stream- 
tube  sufficiently  removed  from  the  heat  release  region  that 
the  flow  is  considered  isentropic.   The  boundary  conditions 
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on  the  bounding  streamtube  are  calculated  by  methods  similar 
to  those  employed  in  determining  the  velocity  and  pressure 
distributions  over  airfoils  of  known  shape.   Since  stagna- 
tion points  do  not  exist  on  the  bounding  streamtube  flow, 
a  modification  to  the  method  presented  by  Sprieter  and 
Alksne  [Ref .  10]  for  airfoils  and  wedges  was  developed. 

In  Section  II,  the  theory  used  in  analyzing  the  transonic 
thermal  blooming  problem  is  presented  and  includes  a  des- 
cription of  the  assumptions  made  during  the  mathematical 
analysis,  a  summary  of  the  pertinent  equations  derived  in 
the  Appendices,  and  a  discussion  of  the  method  of  solution. 
In  Section  III,  the  numerical  results  for  the  Mach  1.0  flow 
with  a  Gaussian  heat  release  distribution  are  presented. 
Lastly,  Section  IV  gives  a  brief  discussion  of  the  conclu- 
sions and  of  the  other  methods  that  show  promise  in  solving 
the  transonic  thermal  blooming  problem. 
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II.   THEORY 

A  laser  beam  that  is  slued  through  the  atmosphere  at  a 
constant  rate  of  ft  radians  per  second  will  experience  regions 
of  subsonic,  transonic,  and  supersonic  flow  in  planes  per- 
pendicular to  its  axis  at  various  radial  positions  along  the 
beam,  Fig.  1.   The  region  to  be  investigated  in  this  paper 
is  the  transonic  region  in  which  mixed  subsonic  and  super- 
sonic flow  exists,  as  in  contrast  to  purely  subsonic  flow  or 
supersonic  flow,  and,  in  particular,  the  case  when  the  Mach 
number  is  precisely  unity;  that  is,  at  a  radial  position 
given  by  Eq.  (1) . 


/7RTc 


rM„  =  i.o  =  r*  =  "IT-  (1) 


For  this  analysis,  the  laser  beam  is  initially  of  a  given 

axial  symmetric  Gaussian  intensity  sluing  into  an  isotropic, 

-3 

quiescent  medium  of  constant  density  p (kgm  m  ) .   Viscosity, 

and  thermal  conduction  are  neglected.   Justification  for 
this  assumption  is  given  in  Appendix  A.   The  heating  effect 
of  the  laser  beam  on  the  medium  is  approximated  by  a  molecular 
relaxation  process  which  is  rapid  enough  that  the  absorbed 
energy  is  instantaneously  transformed  into  heat.   The  energy 
distribution  is  of  the  form  given  in  Eq.  (2) . 

2   2 

Q  =  la  =  IQa  exp(-  (x  +%  }  )  (2) 

2a 
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-3 
where  Q  (Watts  cm   )  is  the  time  rate  of  heat  absorbed  per 

-2 

volume;  Ifl  (Watts  m   )  is  the  beam  maximum  intensity; 

a  (cm   )  is  the  absorption  coefficient;  and  the  exponential 
of  the  Gaussian  with  standard  deviation  a  (cm) ,  which  is 
related  to  the  spot  size  of  the  beam  by  w  =  J 2   a,  is 
centered  at  the  origin  of  the  x-y  coordinate  system  trans- 
verse to  the  beam  axis  at  r*. 

The  theory  employed  for  this  investigation  separates 
the  steady  transonic  flow  with  heat  addition  into  two  en- 
tirely different  flow  regimes,  Fig.  2.   The  internal,  heat 
addition  region  can  be  treated  exactly  by  the  Broadbent 
method  [Ref .  18]  ,  in  which  the  nonlinear  equations  of  motion 
and  the  energy  equation  are  solved  numerically  in  natural 
coordinates.   Then,  a  mesh  of  streamlines  and  normals  is 
constructed  throughout  the  flow  field,  Fig.  3.,  which  reduces 
the  momentum  equations  to  simple  algebraic  equations  solvable 
for  the  flow  properties.   Finally,  the  mesh  of  streamlines 
is  adjusted  until  the  properties  calculated  from  the  linear- 
ized momentum  equations  satisfy  the  energy  equation  and 
boundary  conditions.   Boundary  conditions  are  obtained  from 
the  external,  isentropic  flov;  region  by  approximating  the 
transonic  solution  over  a  bounding  streamtube  to  the  internal 
region. 

The  governing  equations  in  the  internal  region  have  been 
derived  for  flow  with  heat  addition,  but  without  viscosity, 
heat  conduction  or  body  forces.   A  small  perturbation  analysis 
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Cell  boundaries 


FIGURE  3.   Natural  Coordinate  System  with  Flow  Mesh 
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for  transonic  flow  (Appendix  A)  shows  that  viscosity,  heat 
conduction,  and  body  forces  can  be  neglected  to  first  order. 
The  small  perturbation,  steady  two-dimensional  transonic 
equation  with  head  addition  is  written  as: 


(1-M  2U  ,  ,  +d>  ,  ,-M  2(y+l)d>  ,d>  ,  ,  =  77^-   ^{X',y,) 
00  '^x'x1    y  y   °°       yxiyx'x'     CpT   3x' 


(3) 


where  the  term  on  the  right-hand  side  is  added  for  flows 

with  heat  addition;  q(x',y')  (Joules  kgm  )  is  the  heating 

function;  Cp  (Joules  kgm   °K   )  is  the  heat  capacity; 

T   (°K)  is  the  freestream  static  temperature;  M   is  the 
00  00 

freestream  Mach  number;  and  <j>  is  the  nondimensional  velocity 

potential.   Therefore,  in  vector  notation,  the  steady  state 
equations  [Ref.  26]  are: 


V (pU)  =  0  (4) 

pU-vU  +  Vp  =  0  (5) 

pU*V(h  +  jU-U)  =  Q  =  la                       (6) 

p  =  pRT  (7) 


—        -1  -2 

where  U  (m  sec   )  is  the  velocity;  p  (newtons  m  )  is  the 

static  pressure;  h  (Joules  kgm   )  is  the  enthalpy;  and 

R  (Joules  kgm   °K)  is  the  universal  gas  constant.   Using 

a  two-dimensional  natural  coordinate  system  [Ref.  26] , 

shown  in  Fig.  3,  the  governing  equations  were  derived  from 
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the  above  equations.   The  momentum  equations  are  represented 
in  the  following  convenient  numerical  form  (Appendix  B) . 

An' .     .    . . 

ui+l,j+l     ui,j+l  +  ~71 <pi+l,j+l     pi,j+l}         ° 

T  «  (8a) 

YM   2 
Pi+l,j+l-pi,j+l  +    -f-(ui+l,j  +  lvi+l,j+l+ui+l/jVi+l/j)   =0 

(8b) 

where  An!  .  is  the  streamtube  width  at  the  i, jth  station; 
y  is  the  ratio  of  specific  heats;  v!  .is  the  inverse  curva- 
ture defined  as  An,'  ./R;  and  the  prime  (')  denotes  dimension- 

-i- 1  3 

less  quantities  with  respect  to  the  freestream  values.   The 

quantities  with  a  bar  and  parenthesized  subscripts  represent 

average  values  between  the  indicated  stations  (i.e.  An,'.   L,v 

(1,1+D  #3 

is  the  average  streamtube  width  between  station  i,j  and  i+l,j), 
and  the  properties  with  the  tilde  (")    denote  perturbation 
values  from  freestream  conditions  (i.e.  u.  .  =  (u.  .  -  u,  •))  . 

1/3       !  /  3     ■*-  r  J 

The  natural  coordinates  s  and  n  are  tangent  and  normal  to 

the  streamlines,  respectively.   Eqs .  (8a)  and  (8b)  represent 

two  nonlinear  equations  since  u',  p',  An'  and  v1  are  dependent 

variables.   However,  specifying  the  streamlines  of  the  flow, 

these  equations  become  linear  equations  solvable  for 

u'^t  -ji  and  5!.,   . ,  at  point  e  in  terms  of  u!,,  .,  p!,,  • 
1+1,3+1     "1+1,3+1    r  i+l,3   ri+l,3 

at  point  d  and  u I  . , , ,  p !  . . ,  at  point  b ,  Fig .  3 . 
*  i/3+l   ri,D+l    r  3 

The  remaining  flow  properties  p.'  ,  ..,  and  T.'  ,  .,,  can 
3       c      c  ^1+1,3+1       1+1,3+1 

be  calculated  from  the  continuity  equation,  Eq.  (9) ,  and 
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equation  of  state,  Eq.  (10) ,  derived  in  Appendix  B.   How- 
ever, the  new  properties  at  the  point  e  (i+l,j+l)  may  not 
satisfy  the  energy  equation,  Eq.  (11) .   If  it  is  not  satis- 
fied, an  iteration  process  is  performed  by  varying  the 
streamline  distribution  and  recalculating  the  properties 
until  coincidence  exists  throughout  the  entire  flow  field. 
Hence,  the  solution  of  the  problem  is  reduced  to  finding 
the  streamline  distribution  that  satisfies  Eqs .  (8)  through 
(12)  with  the  boundary  conditions  determined  from  the 
external  region. 


pi+l,j+lui+l, j+lni+l,j+l 


=   1  (9) 


Pi+lf  j+1   =    "i+1,  j+lWi+l,  j+1  (10) 

n+i,j+i  -  ;i,j+i +  Ww  (ui+i,j+i  -  ui'j+i>  ■ 

(id 


1    (i,i+l),j    As!/j+1 


Qoo 


where 


p   u  CpT    (l+i(Y-l)M   2)  p   u    (l+i(Y-l)M   2)y 

K0O      00      ^       CO  *  2  00  ^00       00  *  2  '        00  •  .  . 


00 


An,     .  (Y-DAn       . 

-Lf  J  -1-/  J 


Insight  into  the  behavior  of  the  streamlines  at  the 
sonic  speed  was  obtained  from  the  solution  of  the  linearized 
small  perturbation  equation  with  heat  addition  [Refs.  4- 
8] .   A  computer  program  (BLOOM)  was  devised  (Appendix  G) 
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which  gave  solutions  for  a  flow  very  slightly  below  the  sonic 
speed  (Mw  =  0.999),  Figs.  4  and  6,  and  very  slightly  above 
the  sonic  speed  (Mb  =  1.001),  Figs.  5  and  7.   Computation 
of  the  nonlinear  term  in  Eq.  (3)  from  the  linear  results 
indicate  that  the  linear  equations  are  valid  to  M  =  0.9999 
subsonically  and  Mot  =  1.0001  supersonically .   From  the  linear 
solutions,  a  trend  of  the  streamline  shape  is  obtained  in 
the  near  sonic  regime,  even  if  not  at  sonic  conditions. 
Nature  usually  exhibits  a  strong  propensity  to  accomplish 
changes  in  a  rather  smooth  fashion  (even  a  shock  wave  is  a 
smooth  change  in  conditions  if  viewed  microscopically 
enough);  thus,  a  knowledge  of  flow  conditions  in  the  near 
sonic  range  gives  a  good  clue  to  the  sonic  flow  behavior. 
In  Figs.  4  through  7,  the  Gaussian  heat  distribution  is 
centered  at  the  origin  on  the  centerline,  and  the  relative 
airflow  is  from  left  to  right. 

From  the  appropriate  density  calculations  and  plots, 
discussed  later  in  the  thesis,  it  can  be  seen  that  the  wake 
term  becomes  less  significant  compared  with  the  source  term, 
which  is  directly  associated  with  the  horizontal  velocity 
perturbation,  as  sonic  velocity  is  approached  from  either 
direction  (i.e.  M  <  1  and  M  >  1) .   Figs.  4  and  5  show  the 
horizontal  (longitudinal)  velocity  perturbation  at  Mach 
numbers  of  0.999  and  1.001,  respectively.   Figs.  6  and  7 
show  the  vertical  (transverse)  velocity  perturbation  at  these 
same  Mach  numbers.   Note  that,  subsonically,  the  horizontal 
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perturbation,  Fig.  4,  decreases  to  a  negative  valley 
slightly  upstream  of  the  laser  beam  center  and  then  changes 
rapidly  to  an  equal  positive  peak  slightly  farther  downstream 
of  the  beam  center.   The  perturbation  contours  are  of  very 
small  curvature  near  the  beam  center;  the  supersonic  per- 
turbations, Fig.  5,  show  a  negative  valley  near  the  beam 
center  with  a  return  to  freestream  conditions  downstream  of 
the  center.   In  the  supersonic  case,  the  contours  are  nearly 
normal  to  the  flow  throughout,  in  keeping  with  the  character 
of  supersonic  flow;  whereas,  in  the  subsonic  case,  the  con- 
tours are  closed,  in  keeping  with  subsonic  flow  behavior. 
The  patterns  of  disturbance,  however,  appear  to  be  consistent 
as  the  sonic  speed  is  crossed. 

The  vertical  or  transverse  velocity  perturbation,  Figs. 
6  and  7,  show  similar  behavior  in  the  two  speed  regimes. 
As  expected,  the  subsonic  flow  perturbation  profiles  are 
closed  while  those  for  supersonic  flow  are  open  to  infinity. 
In  both  cases,  the  perturbations  exhibit  an  odd  symmetry 
about  the  main  flow  centerline.   For  the  subsonic  case,  the 
perturbations  are  greater  than  zero  to  the  left  of  the  rela- 
tive flow  and  less  than  zero  to  the  right.   This  is  a  direct 
indication  of  the  spreading  of  the  streamlines  as  the  heat 
release  region  is  approached  by  the  air  flow.   Downstream 
of  the  laser  beam,  the  streamlines  again  become  parallel 
with  the  undisturbed  flow  direction.   The  supersonic  case 
shows  the  same  patterns  except  that  the  spreading  effect  on 
the  streamlines  is  felt  infinitely  far  in  the  transverse 


31 


direction,  roughly  along  the  characteristic  lines  from  the 
disturbances. 

The  effect  of  spreading  of  the  streamlines  as  the  beam 
is  passed,  which  can  be  constructed  from  these  plots  of 
longitudinal  and  transverse  velocity  perturbations,  is  thus 
seen  to  be  quite  similar  immediately  below  and  immediately 
above  the  sonic  speed.   Fig.  2  is  a  sketch  of  the  stream- 
lines.  It  is  expected  therefore  that  the  streamline  behavior 
at  the  sonic  speed  will  not  vary  radically  from  this  pattern. 

Based  on  the  streamline  distribution  in  Fig.  2,  the 
following  relationship  between  area  change  (streamtube  width 
variation)  and  change  in  total  temperature  (heat  input)  was 
chosen,  Eq.  (13) . 


dA   dn    1  #.,*/,.«, »   o 


x  -  t  -  f(Y+1)  (1+s)^ 


(13) 


where  T  is  the  local  total  temperature  of  the  flow  and  3 
is  a  variable.   One  dimensional  influence  coefficients  with 
area  change  and  heat  addition  [Ref .  27] ,  Eq.  (14)  ,  were 
used  to  determine  the  required  $  at  precisely  Mach  1.0. 


dM2  =   G(x) 

dx  '   (1-M2) 


(14) 


where 


?        i  9  HirWr^  9    dln(T    ) 

G(x)    =   VT(l+±{y-l)vf)  i-2    ax  +  <1+?M   >     dx  } 
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In  order  for  Eq.  (14)  to  remain  finite  and  to  be  defined 
at  Mach  1.0,  G(x)  must  tend  to  zero  in  the  limit  as  the  Mach 
number  approaches  unity.   If  Eq.  (13)  is  substituted  into 
Eq.  (14)  and  the  limit  taken,  it  is  found  that  it  must  be 
zero  at  precisely  Mach  1.0  (Appendix  C) .   Hence,  the  stream- 
tube  area  change  is  related  to  the  heat  input  by  Eq.  (15) . 


dn!  . 
1/3 


dT" 
o 


-£<Y+1>   _id 


'1,3 


(15) 


where  the  maximum  difference  between  n!    and  n'  .,  and  T' 

i,:       1,3  o.  . 

-6  'J 

and  T"     is  of  the  order  of  10    and,  hence,  the  freestream 
o,  . 

values  can  be  used  in  the  denominator.   Using  this  approxi- 
mation for  the  streamtubes,  an  explicit  formulation  of  the 

streamtube  shape,  at  any  position  s!  .  along  a  streamtube, 

1 1  3 

was  obtained  (Appendix  C) . 


n 


;,i  =  i(AHi,i+1);   i^1 


(16a) 


n 


i/D 


1        ^ 
=  ±An!  ,  +   Z   An! 


2   i,l 


1=1 


i,l+l 


+  (J  " 


§><• 


i  >  1,  j    2 


(16b) 


where  for  a  Gaussian  heat  input  at  Mach  1.0  (Appendix  C) , 


An!  . 
i/D 


20 exp("   .2   } 

^w°o  2a' 


s! 

1+erf  ( — — ) 
72a' 


;  i   1/  j   1 


(16c) 
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Once  the  streamtube  shape  is  known,  the  slope  (dy/dx) , 
Eq.  (17) ,  and  the  radius  of  curvature  (R) ,  Eq.  (18) ,  can 
be  calculated  at  every  mesh  point  throughout  the  flow. 


£      =  6      =  ni+i.rWti,.i-ni,i+frni,i.  i  >  x 


i  >  i 

(17) 


*»*      R     AS!  .(l+«2  .)3/2        isi,j      '    -' 

(18) 


Therefore,  the  proportionality  between  the  streamtube 
shape  and  heat  release  distribution  gives  not  only  an 
explicit  formulation  for  the  streamtubes  and,  hence,  linear- 
ization of  the  nonlinear  momentum  equations  but  also  simul- 
taneous agreement  between  those  properties  calculated  from 
the  momentum  equations  and  the  energy  equation.   The  problem, 
in  the  case  of  precisely  M  =1.0  flow,  becomes  one  of 
determining  the  correct  boundary  conditions. 

The  boundary  conditions  that  are  used  are  those  corres- 
ponding to  freestream  conditions  (ul     .  =  p'  .  =  p '  .  =  T,*  •  -  0)  , 

-L/J       -"-fJ       -LfJ       -wJ 

which  are  far  enough  upstream  that  the  properties  are  not 
affected  by  the  perturbation  caused  by  the  heat  release 
region,  and  the  pressure  (p!  .  ic)    and/or  velocity  (u!  .*)  on 
the  j*th  streamtube  which  is  considered  far  enough  removed 
from  the  heat  release  region  that  the  flow  can  be  considered 
isentropic.   Eqs.  (16b)  and  (16c)  can  be  used  to  determine 
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the  streamtube  shape  and  displacement  where  j  =  j*.   However, 
a  more  convenient  form  is  derived  in  Appendix  C  for  a 
Gaussian  heat  release  distribution  at  M  =1.0  which  can  be 

00 

written  as: 


n! 


i<D 


-n' 


2 

(y  +  1)  I  ao  '  7T 


l/j 


4Q 


s! 
1+erf  ( — — ) 


J2a 


;  i  >_  1,  j  >_  4.5 
(19) 


and  the  slope  of  the  bounding  streamtube: 


dn!  .* 

*s:,j*  ■ 


(Y+l)I0a  J2^a% 


,2 


4Q 


-exp  - 


2a 


7? 


f        f 
-  exp 

2iro '      2a 


U 


(20) 


where  t1  is  the  maximum  growth  in  the  streamtube. 

Since  the  mass  flow  is  constant  in  a  streamtube,  the 
bounding  streamtube  can  be  considered  as  a  solid  wall  simi- 
lar to  the  surface  of  a  wedge  with  a  cusp  extending  forward 
to  the  freestream  conditions  upstream  at  infinity,  Fig.  2. 
By  hypothesis,  this  flow  is  isentropic  and  can  be  treated 
by  methods  similar  to  those  used  in  analyzing  transonic 
airfoil  and  wedge  flows.   In  particular,  an  approximate 
solution  is  obtained  through  an  iteration  process  on  the 
integral  equation  derived  by  Green's  functions  for  small 
disturbance  transonic  flow  similar  to  that  performed  by 
Sprieter  and  Alksne  [Ref.  10  and  11]  for  thin  airfoils  and 
wedges  but  modified  for  the  streamtube  configuration  which 
has  no  stagnation  points .   The  absence  of  stagnation  points 
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in  the  case  of  the  flow  around  the  bounding  streamtube 
enhances  the  accuracy  of  the  approximate  numerical  tech- 
niques because  the  stagnation  point  represents  a  singularity 
in  the  problem  and  places  added  constraints  on  the  method  and 
accuracy  of  the  solution  in  the  region  of  the  stagnation 
point.   In  this  method,  the  quadratic,  nonlinear  nature  of 
the  governing  equations  is  retained  which  precludes  shock- 
free  supercritical  flows  in  the  transonic  regime. 

The  integral  equation  for  transonic  flow  is  derived  by 
a  Green's  function  analysis  on  the  small  perturbation 
transonic  equation  without  heat  addition,  Eq.  (3)  without 
the  right-hand  side,  resulting  in  Eq.  (21) . 


uw(x,0)  =  uL  (x,0)  +  |u^(x,0)  -  ~I  (21) 

w 


where 


M   (y+1)u*  (x,0)    M  -  M 

u  (x,0)  =  — =  5_^_  (22) 

6^  3 


uLU,o)  -i   /  **_£_  (23) 

W  -oo      ^  (x-£) 

2  — 
oo  u  (x,0)    —    — 

1=    /   -E-5 E(*  g  X)  d£,  (24) 


—  00 


—   —       -    oo   2     r, 
E(^  7  X)  «  -  £-      /   2_  in  -I  dn  (24) 

b         27T  0    dC  r2 
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and  the  normalized  variables  are  defined  in  Appendix  D. 
Eq.  (21)  is  a  function  only  of  x,  the  normalized  coordinate 
along  the  wall  surface,  and  t,  the  normalized  amplitude  of 
the  wall  slope,  given  by: 


where 


—  —2 

dZ     —  X                                            /oe  \ 

—  «  t  exp  -  —^                                                                       (25a) 

dx  2o* 


Mm2  (y+1)  t   Mm2  (y  -1)  la  J2T  a 

t  =  . =  . 2 (25b) 

Z  S  P  u  y 


The  iteration  technique  (Appendix  D)  is  performed  for 
various  subsonic  freestream  Mach  numbers  starting  at  the 
initial  supercritical  flow  and  proceeding  toward  the  Mach 
1.0  flow.   Each  solution  obtained  by  the  iteration  technique 
gives  a  unique  T,  which  determines  the  freestream  Mach  number 

and  its  corresponding  normalized  velocity  distribution 

—  —  —  —2/3 

u  (x,0)  along  the  surface.   If  (u  -1)  versus  t    is  plotted 

for  each  position  along  the  surface,  it  is  found  that  as 

the  freestream  Mach  number  tends  toward  unity  the  slope 

becomes  constant.   This  corresponds  to  the  phenomenon  of  the 

Mach  number  freeze  [Ref.  2  6]  wherein  the  local  Mach  number 

is  invariant  with  changes  in  the  freestream  Mach  number 

when  the  latter  is  near  unity  or,  more  precisely, 


ar        =  °  <26' 

00  M  =1 
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Vincenti  and  Wagoner  [Ref.  28],  Liepmann  and  Bryson  [Refs. 
29  and  30] ,  and  others  have  shown  that  the  corresponding 
approximate  relation  yielded  by  the  small-disturbance 
transonic  theory  is  given  by: 


3i 

at 


=  0  (27) 


where 


E (1-M2)      _  (u-1)  r2p, 

C  "    /u  2,  ~  ,2/3 =2/3  (28) 

M   (y+1) t)  t 


oo 


The  freeze  of  Mach  number  extends  over  a  finite  range  near 
Mach  1.0  where  it  is  constant.   Hence,  £,  the  slope  of  the 
(u-1)  versus  t  curve  at  each  location  x,  is  constant  near 
the  sonic  condition,  and  the  local  Mach  number  at  each 
position  can  be  calculated  by  Eq.  (28) .   With  the  local 
Mach  number  distribution  known,  the  velocity  perturbation 
can  be  calculated  from  Eq.  (22) ,  and  the  other  flow  proper- 
ties can  be  determined  from  isentropic  flow  relations  [Ref 
27]. 

The  Mach  1.0  flow  is  completely  specified.   Boundary 
conditions  are  known  upstream  of  the  heat  release  region 
and  on  a  bounding  streamtube.   Calculation  of  the  entire 
flow  field  is  readily  calculated  from  the  simple  algebraic 
equations  obtained  by  Broadbent's  method,  where  the  stream- 
tube  configuration  is  known  precisely  at  Mach  1.0. 
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III.   NUMERICAL  RESULTS 

The  computed  flow  properties  are  for  the  Mach  1.0  thermal 
blooming  problem  with  the  characteristics  listed  in  Table  I. 


Gaussian  Beam 

Peak  Intensity   (I  ) 

o 

Standard  Deviation   (a) 

3.333xl07  Watts /m2 
5  cm. 

Atmospheric  Absorption   (a) 

-7    -1 
8.0  x  10    cm 

Atmospheric  Freestream 
Conditions 

Temperature   (T  ) 

Pressure   (p  ) 

Density   (pj 

Velocity   (u  =  a*) 

Streamtube  Width 
<Anl,j=AnJ 

Streamline  Cell  Length 
(AS.(1) 

288.0   °K 

1. 013x10 5  Newtons/m2 
1.2246   kgm/m 
340.0  m/sec 

1  cm. 

1   cm. 

Ratio  of  Specific  Heats   (y) 

1.4 

Specific  Heat  at  Constant 
Pressure   (Cp) 

1005.0   Joules/kgm°K 

Table  I 

Freestream  Flow  Properties  And 
Laser  Beam  Characteristics 
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The  first  step  in  the  procedure  is  to  determine  the 
boundary  conditions  for  the  heat  release  region.   Far  up- 
stream of  the  heat  release  region,  all  the  dimensionless 
perturbation  quantities  are  zero.   The  other  boundary  con- 
ditions are  calculated  by  Sprieter's  method  [Ref.  10]  for 
the  bounding  streamtube,  the  shape  of  which  depends  upon 
the  heat  release  distribution  and  is  known,  Eq.  (16) . 

Using  Sprieter's  numerical  iteration  procedure  on  the 
transonic  flow  integral  equation,  Eq.  (21) ,  the  normalized 
velocity  distribution  along  the  bounding  streamtube  was 
obtained  for  several  frees tream  Mach  numbers  between  the 
initial  supercritical  flow  and  Mach  1.0.   The  results  are 

graphed  in  Fig.  8.   From  these  solutions,  a  (u-1)  versus 

—2/3 

t   '  curve,  Fig.  9,  was  plotted  for  several  positions  along 

the  streamtube.   The  results  indicate  that  the  slope  (£) 

is  constant  for  each  position  and,  hence,  the  Mach  number 

can  be  considered  frozen. 

From  the  values  of  Clx),  the  local  Mach  number  variation 

at  a  freestream  Mach  number  of  unity  was  calculated  from 

Eq.  (28)  ,  since  %   is  known  and  constant  for  the  known 

physical  streamtube  shape,  Eq.  (30) ,  from  Eq.  (D6) . 


x  =  — - —  =  1.386  x  10"6  (30) 

J2T?  a 


where 


2 
(Y-l)Lcnra     ,  -,->-,    ,  A-5 
t  =  — ■ =  1.737  x  10    cm. 


YP  u 

1  *■  OO   00 
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(a)  Shock  Wave  at  x  =  9.5 


-10     l  -5 
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1/e  of  Gaussian  beam 


u 


(b)   Shock  Wave  at  x  =  11.5        j.,5  • 
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(c)   Shock  Wave  at  x  =  14.5 


(e)   Shock  Wave  at  x  =  24.5 
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FIGURE    8.      High  Subsonic   Supercritical   Flows   on 
Bounding  Streamtube 


41 


1.0 
0.8 

0.6 

0.4 

0.2 

(u-1)  0 

-0.2 

-0.4 

-0.6 

-0.8 

-1.0 

-1.2 

-1.4 

-1.6 

-1.8 
-2.0 


(10.5) 


(3.5) 


(2.5) 


(1.5) 


(0.5) 

(-29.5) 
(-24.5) 

(-19.5) 


(-9.5) 
(-4.5) 


*Note:  Number  in  parentheses  represents  x  position.         ' 


FIGURE   9.      Mach  Number  Freeze   £   =   Constant   for 
Various    x  Positions 
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is  the  maximum  growth  of  the  bounding  streamtube. 

Once  the  local  Mach  number  was  determined,  the  velocity 
perturbation  along  the  bounding  streamtube  was  calculated 
from  Eq.  (22)  and  is  plotted  in  Fig.  10.   The  remaining 
flow  properties  were  calculated  by  the  isentropic  flow 
relations  [Ref .  27] . 

It  has  been  shown  that  the  area  change  of  the  stream- 
tubes  in  the  heat  addition  region  is  directly  proportional 
to  the  change  in  total  temperature;  that  is,  6  =  0  in 
Eq.  (13) .   The  energy  equation  is  automatically  satisfied 
by  this  selection  for  the  area  variation.   Since  the  heat 
release  distribution  is  known,  all  of  the  streamtubes  are 
specified.   From  the  streamtube  shape,  the  radius  of  curva- 
ture of  the  streamlines  were  calculated  at  every  point 
throughout  the  region,  Eq.  (18) . 

With  the  curvature  known,  the  calculation  of  the  flow 
properties  in  the  heat  region  is  reduced  to  the  solution 
of  two  simple  linear  algebraic  equations,  Eqs.  (8a)  and 
(8b),  the  momentum  equations  in  natural  coordinates.   The 
results  of  this  "marching"  technique  are  shown  in  Fig.  11 
where  the  density  perturbations  in  the  upper  half  of  the 
symmetric  laser  beam  are  presented.   This  method  has  the 
distinct  advantage  that  it  can  be  solved  on  an  HP  9830 
computer . 

Figs.  12  and  13  show  the  linear  supersonic  and  subsonic 
density  perturbation  solutions  calculated  from  the  computer 
program  (BLOOM)  developed  in  Appendix  G.   The  results 
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presented  in  Fig.  10  at  Mach  1.0  are  consistent  with  a 
continuous  transition  through  the  transonic  regime.   Although 
the  solution  calculated  is  for  precisely  Mach  1.0,  it  can 
be  logically  hypothesized  that  the  flow  behavior  through 
the  transonic  regime  is  as  shown  graphically  in  Fig.  14. 

The  linear  supersonic  and  subsonic  solutions,  and 
Sprieter's  high  subsonic  solutions  are  in  excellent  agree- 
ment with  this  hypothesis  as  well  as  that  conjectured  from 
the  hodograph  plane  representation  for  transonic  flow 
(Appendix  E) .   Therefore,  there  is  a  steady,  two-dimensional 
solution  for  the  transonic  thermal  blooming  problem  which 
is  finite  and  does  not  result  in  "catastrophic"  defocusing 
of  the  laser  beam. 
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IV.   CONCLUSIONS 

It  has  been  demonstrated  that  there  is  a  two-dimensional 
steady  state  solution  to  the  transonic  thermal  blooming 
problem  for  a  sluing  laser  beam  with  a  freestream  Mach  number 
of  precisely  unity.   According  to  the  linearized  solutions 
for  thermal  blooming,  the  density  perturbations  become 
infinite  as  a  Mach  number  of  unity  is  approached.   Due  to 
nonlinear  effects  the  trend  to  infinity  is  cutoff,  and  the 
finite  value  of  perturbation  quantities  has  been  established. 

In  obtaining  the  solution,  it  was  found  that  at  precisely 
Mach  1.0  that  the  streamtube  area  variation  through  the  heat 
release  region  is  directly  proportional  to  the  change  in 
total  temperature  in  that  region.   With  the  streamtube 
shapes  known,  an  exact  solution  for  the  flow  field  can  be 
obtained  by  Broadbent's  method  [Ref.  21]  with  the  corres- 
ponding boundary  values  for  Mach  1.0  isentropic  flow  far 
from  the  heat  release  region.   There  are  no  restrictions  on 
the  existence  of  shock  waves  nor  on  the  symmetry  or  asymmetry 
of  the  heat  release  distribution,  provided  the  bounding 
streamtube  in  the  isentropic  region  and  the  boundary 
conditions  on  it  can  be  calculated. 

Furthermore,  in  obtaining  the  solution  at  precisely 
Mach  1.0,  it  was  found  that  shock  waves  exist  throughout 
the  transonic  regime,  Fig.  14.   As  the  sonic  condition  is 
approached  from  subsonic  speeds,  a  shock  wave  forms  at  the 
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downstream  1/e  (one  spot  size)  size,  of  the  Gaussian  heat 
release  distribution,  where  the  local  Mach  number  first 
becomes  sonic,  and  rapidly  moves  downstream  to  infinity. 
If  the  sonic  speed  is  approached  supersonically,  a  shock 
is  formed  at  the  beam  center.   To  be  consistent  with  the 
results  in  the  subsonic  portion  of  the  transonic  regime, 
it  has  been  hypothesized  that  the  shock  moves  rapidly  up- 
stream to  infinity.   Finally,  the  analysis  in  Appendix  F 
indicates  that  there  is  a  minimum  laser  sluing  rate  that 
must  be  maintained  to  preclude  significant  phase  distortion 
in  the  beam  at  Mach  1.0. 

The  results  of  the  Mach  1.0  thermal  blooming  problem 
for  a  Gaussian  intensity  distribution  are  in  excellent 
agreement  with  the  approximate  results  obtained  by  Ellinwood 
and  Mirels  [Ref.  19].   Fig.  15a  shows  these  results  where 

the  maximum  density  perturbation  was  calculated  to  be  of 

2/3 
the  order  of  Q  '  .   Q  is  defined  as: 


(y-1)  I  ar 
r 

ypa 


Q  =  —  ,  (31) 


where  r  is  the  radius  (spot  size)  of  a  circular  beam  of 
uniform  heating  intensity  (I  o°0  •   This  approximation,  when 

applied  to  the  Gaussian  beam  of  Section  III,  gives  a  maximum 

-4 
density  perturbation  of  the  order  of  10   which  is  consistent 

with  that  obtained  in  this  paper,  and  shown  in  Fig.  15b. 

Although  the  density  perturbations  are  very  small  (i.e. 

measured  in  parts  per  million) ,  significant  optical  defocusing 


51 


^2/3 


'^-Linear  Theory 

/  \ 
/  \ 


Nonlinear  Theory 


Ap 
P 


0  1.0 

Mach  Number,  M  =  far/a^ 
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(b)   Results  of  Linear  Subsonic  and  Supersonic  Solutions  from  DIDER 
[Refs.  4  and  5]   and  of  Broadbents  Method  at  Mach  1.0 

FIGURE   15.      Schematic   Representation  of   the   Steady   State 
Density   Perturbations    for   a  Sluing 
Two-Dimensional   Laser   Beam 
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of  the  laser  beam  could  result  (Appendix  F) .   With  an 

increase  in  the  beam  area  due  to  defocusing,  a  large  decrease 

2 

in  the  beam  intensity  I  (Watts/m  )  will  occur. 

Two  other  methods  that  show  promise  in  solving  the  tran- 
sonic thermal  blooming  problem  are  the  finite  difference 
method  developed  by  Murman  and  Cole  [Ref.  12]  for  isentropic 
transonic  flow  which  can  be  extended  to  include  heat  addi- 
tion, and  the  "finite  element  method"  of  solving  partial 
differential  equations  [Ref.  31]  which  is  a  powerful  method 
for  solving  a  wide  range  of  engineering  problems . 

A  Mach  1.0,  steady  state  solution  for  thermal  blooming 
in  a  laser  beam  with  a  given  Gaussian  heat  release  distribution 
has  been  presented  in  this  thesis,  Fig.  11.   Subsequent  inves- 
tigations should  be  performed  to  ascertain  the  Mach  1.0 
solutions  for  various  Gaussian  and  non-Gaussian  heat  release 
distributions,  beam  intensities  (I),  atmospheric  absorption 
coefficients  (a),  and  frees tream  flow  properties.   In  this 
manner,  a  correlation  might  be  developed  that  would  be  use- 
ful as  an  engineering  approximation  to  the  density  perturba- 
tions at  various  transonic  Mach  numbers,  vice  a  simple  order 
of  magnitude  linear  extrapolation  as  presented  by  Ellinwood 
and  Mirels  [Ref.  19] .   Although  solutions  are  known  at  Mach 
1.0  and  in  the  linear  subsonic  and  supersonic  regions,  from 
the  computer  program  BLOOM,  it  would  be  premature  to  estab- 
lish a  criterion  that  could  be  used  in  estimating  transonic 
behavior  in  general  before  additional  work  is  completed. 
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APPENDIX  A 

DERIVATION  OF  STEADY  TRANSONIC  FLOW 
FOR  A  CONDUCTING  GAS  WITH  HEAT  ADDITION 


Eqs.  (Al)  through  (A4)  govern  steady  compressible  flow 
with  heat  addition  and  constant  properties  [Ref.  32],   These 
equations  include  viscosity  and  heat  conduction.   The 
distribution  of  the  heat  addition  is  given  by  the  function 
Q  f(x,y)  and  is  assumed  known.   The  gas  flow  is  assumed 
uniform  at  infinity  (x  =  ±«) . 

V (pU)  =  0  (Al) 

pU-VU  +  Vp  =  ip  (A2) 

where       ty  e  nV(V«U)/3  +  nV2U 

pTU«Vs  ■  -  fc-  Vq  (A3) 

where       <D  =  —p-  (V  -U)  2  -  n(V2(U-U)  -  (VXU)2  -  2U-V2U) 

2 
and         Vq  =-kV  T  -  Q  f(x,y) 


2  2 

Vp  =  pa  Vs/Cp  +  a  Vp  (A4a) 


or 


s  -  s 
expt-^^-0-)  =  T/p(Y_1)  (A4b) 
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The  following  dimensionless  variables  are  defined: 


x'  =  x/l1 

y'  =  y/^2 

P'   =  P/Pco 

T1  =  T/T 

'    00 

U'  =  U/U^ 

S«  =  s/Cv^ 

n'  =  n/nB  =  1 

K '  =  K/K   =  1 

'    00 

ere  £,  and  l~   are 

characteristic 

lengths  in  the 

directions,  respectively,  and  the  reference  values  are  the 
freestream  conditions.   Substituting  the  dimensionless 
variables  into  Eqs .  (Al)  through  (A4)  and  eliminating  the 
pressure  term  in  the  momentum  equation,  Eq.  (A2) ,  with  the 
equation  of -state,  Eq.  (A4) ,  the  equations  simplify  to  the 
following: 

p» V •U'    +   U,-7,-pl    =   0  (A5) 


p.P.. 7-o»  =  -  sjl§-  - p  s-y  *    ;vRe°   +  v-    (A6) 

Moo  ^Moo 

Y(Y-DMM2    (2  _      ;  ,   _     _ 

p'T'U'-V's'    =   — <    -   j    (V'U1)^    +    V,Z(U'-U) 


-    (V'XU1)2    -    2(U' •V,2U,)j 


+   YRepr'    +   Y(Y-DMw2Haf(x,y)  (A7) 


where  the  dimensionless  ratios  are  defined  as  follows 
Reynolds  number:      Re  =  p  U  1-,/r) 


Mach  number:  M  =  U  /VyRT 


00        oo 


Prandtl  number:      Pr  =  C   n  /< 

n   'oo'   oo 

00       3 
Heat  quantity:        H   =  Q  £,/p  U 

-1  J  00        O   X     00   00 

Scaling  factor:      L  =  l^/l~ 
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For  flows  in  which  conduction  and  viscosity  can  be 
neglected,  two  flows  will  be  dynamically  similar  when  y, 
M   (Mach  number  in  the  absence  of  heat  addition)  and  H 

00  00 

are  the  same. 

Retaining  the  viscous,  conduction,  and  heat  addition 
terms,  the  governing  equations,  Eqs .  (A5)  through  (A7),  in 
two-dimensions  become: 


,  ,3u'    ,  3v\     ,  3p'    _  ,  ipj      n 

p  (33T  +  L  W]  +  u  &  +  Lv  ay1"  =  ° 


(A8) 


x-direction  momentum  equation: 


,  ,  ,  3u'  ,  _  ,  3u ' x 
p1  (u1  -r— r  +  Lv'  -r--r) 
3x'       3y* 


2„, 


T'  3p*  _  p'T'  8s'     4   3  u 


M 


2  3x 


YMr 


2  3x 


3Re  ~  ,2 
3x' 


3  v1 


L2  32u' 


+  3Re  3x'3y'  +  Re  "^~T2~ 


(A9) 


y-direction  momentum  equation: 


i  /    i    3V        _     ,    3v' « 

p    (u      3F-  +   LV      W] 


LT'     3p'         Lp'T'     3s'         4L2    32y' 

P1^'   YM  2  3y     3Re^ 


_L_      32u'      +     1     32v' 


3Re    3x'3y'         Re    ~    ,2 


(A10) 


energy   equation: 


p'T-tu-  h;+lv  aa;, 


•        ■  c 

Re 


4    f  ,3u\  2  3u'3v' 

3   1 ^3x';  ^   3x'3y' 


32V 


3y 


,2 


,     f3v«         _    3u'2 


^    (3x'2 


,2 


Y(y-l)MmzHm   f(x,y) 


(All) 
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Since  the  transonic  flow  region  is  of  interest/  a 
transonic  expansion  procedure  similar  to  that  used  for 
incompressible  flow  around  an  airfoil  of  thickness  ratio  t 
[Ref.  33]  will  be  used.   In  the  following  analysis,  the 
heat  addition  will  cause  the  perturbations  from  the  free- 
stream  conditions.   Introducing  the  following  perturbation 
quantities : 

u'  =  1  +  e2/3  u^  +  e4/3  uj;   ;  v'  =  ev[ 

p'    =   1   +    e2/3    p»    +    e4/3    p£      ;  T'    =   1   +    e2/3   T{   +    e4/3   T£ 

s'    =   1  +   e4/3   s£      ;           1/Re    -  0(e)       ;           Hoo~0(e4/3)       ; 

lx  =  1        ;           l2  =   s"1/3        ;  L  =   s1'3 

the  first  and  second  order  equations  at  M^  =  1.0  become: 

First  order  equations: 

Su,'    8p\' 

j£   +  8^  =  °  (A12) 


3v'    3p' 

1ST  +  3^  "  °  (A13) 


3T'  8u' 

__1  +  (Y-D  __1  =  o  (A14) 


Integrating  Eq.  (A12)  gives  ui  =  -pi    which,  when  substituted 
into  Eq.  (A13) ,  illustrates  that  for  first  order  theory  the 
flow  is  irrotational. 
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3v'    3u' 

I3F  -  W  m   °  (A15> 


Finally,  Eq.  (A14)  gives  the  temperature  distribution  as: 


TJ  =  (Y-D  P{  =  -(Y-D  u[  (A16) 


Second  order  equations: 


3ul  3u'         3v'         3p'  3p' 

1HF  +  pi  a^  +  3y^  +  a^  +  ui  33F  -  °  (A17) 


3u£  3u^  3u.[  x        3p£  3p| 

33Tr+ui3Tr+pi3lTr="TT  (3P-  +  Tl  JST) 

M 

00 

.       3SI        4e1/3   32u' 

-  -rr  3^  +  -r-  r^         (A18) 

YM  3x' 

'        00 

3vJ  3vJ  x        _      3p{        3p£ 

pi  3P"  +  ui  33F  =  "  ^T  (Ti  371-  +  ay1"1 


CO 


1      3§2   +      1/3    ^i 


YM   2    3y'  3x'2 


(A19) 


8§2  1/3    3    *1  2 

__4  =  ye-L/J  1  +   Y(Y_i)    M^   f (x,y)  (A20) 

3x' 


Combining  Eqs .     (A17)    through    (A19)    and  using  Eqs .     (A12)    and 
(A16) ,    the   transonic   perturbation  equation   is   obtained. 

3u'         3v»  .    1/3    32u'  ,  /r,    32T' 

-(Y+l)    „•  ^  +  ^i  -  -  ifi,-        1  +  EV3         1  +   (y.1)f  (x#y) 

d  X  o  X 

(A21) 
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The  first  term  on  the  right-hand  side  of  Eq.  (A21) 
represents  the  effect  of  viscosity,  the  second  term  the 

conduction,  and  the  last  the  heat  addition.   As  can  be  seen, 

1/3 
both  the  viscous  and  conduction  terms  contain  a  e  '  "  term 

when  the  Reynolds  number  is  of  the  order  of  e   ,  and,  hence, 

are  small  and  can  be  neglected  under  most  circumstances. 

-2/3 

If,  however,  the  Reynolds  number  is  of  the  order  of  e 

or  smaller,  both  terms  become  important  and  must  be  retained 
in  the  perturbation  equation.   Since  to  first  order  the  flow 
is  irrotational,  the  velocity  potentials  u|  =  34>/3x'  and 
vl    =  9<j>/3y'  can  be  introduced  into  Eq.  (A21)  giving: 

(3y+1) e1' 3 
-CY+1)  V»x'x'  +  Vy'  ="  i2LH Vx'x' 

+  (Y-D  f  (x',e"1/3y')        (A22) 


A  similar  analysis  can  be  carried  out  for  a  freestream 

Mach  number  close  to  unity.   In  this  case,  the  freestream 

2/3 

Mach  number  can  be  approximated  by  Ma  =  1  +  Ke    and  the 

transonic  expansion  gives: 


1/3 
(K+(Y+l)*xJ*x'x«  "Vy  =  (3Y+3)£     *x'x'x' 


-  (Y-l)  fU'^'^V)   .   (A23) 
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APPENDIX  B 

BRIEF  DEVELOPMENT  OF  BROADBENT'S  METHOD  OF 
SOLVING  NONLINEAR  EQUATIONS'  OF  MOTION  WITH 

HEAT  ADDITION 


Broadbent's  method  [Ref.  24]  is  an  inverse  solution  in 
which  the  streamlines,  normals,  and  boundary  conditions 
along  at  least  one  streamline  and  normal  are  specified  and 
the  resulting  flow  field  calculated. 

The  governing  equations  for  this  derivation  are  as 
follows : 

V(pU)  =  0  (Bl) 

pU*VU  +  Vp  =  0  (B2) 

pTJ.V(h  +  jU*U)  =  Q  =  la                       (B3) 

p  =  pRT  (B4) 

In  natural  coordinates  [Ref.  26] ,  these  equations 
become : 


|-  (puAn)  =  0  (B5) 

dS 


pu  If  ♦  If  =  0  (B6) 

P  4   +  f  =  °  (B7) 

pu  |^  (h  +  iu2)  =  la  (B8) 


60 


A  mesh  is  constructed  over  the  flow  field  using  stream- 
lines and  normals  which  are  chosen  in  advance  as  in  Fig.  3. 
The  cells  specify  the  streamtube  boundaries  and  the  mesh 
lines  the  streamline  slope  and  curvature  throughout  the 
flow  field.   With  the  mesh  specified,  Eqs .  (B5)  through  (B7) 
are  nondimensionalized  using  uniform  freestream  flow 
properties  (u,  . ,  p,  . ,  p,  . ,  and  T..  . )  ,  a  constant  uniform 
streamtube  width  upstream  of  the  perturbations  (An,1  .)  and 

■*-  r  J 

constant  Cp. 


o!  .u!  -An!  ,  =  1  (B9) 

Pi, 3  i,3a  1,3 

u'       -u!      =  -  An^-i+1)^  (p.       -p-     ) 
ui+l,j+l   Ui,j+1        yM  2      ^pi+l,j+l  Pi,j+lj 

(BlOa) 


*i+l,j+l-Pi+l,j  =  "  I  YM«2(ui+l,j+lVi+l,j+l 


+  u'    .v].,  •)        (BlOb) 
1+1,3  1+1/3 


n+i,j+i-n,j+i  ^"'^""w.h-O  = 


q'  '   ",J  isi,j  +  l  =  OI.j+1    (B11) 


p!  .  =  p!  .  RT!  . 
*i,3    Mi,3    1,3 


(B12) 


where  Qw  is  a  freestream  heat  quantity  defined  in  Eq.  (12) ; 
Qf   ,  .  is  a  dimensionless  heating  term;  and  v.'  .  is  a 
dimensionless  inverse  radius  of  curvature. 
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The  slope  (dy/dx)  and  the  radius  of  curvature  (R)  of 
the  streamlines  at  a  point  b  (i+l,j)  are  shown  geometrically 
in  Figs.  3.   Due  to  the  symmetry  of  the  problem,  the  slope 
and  the  curvature  of  the  centerline  streamline  are  zero 
(i.e.  for  all  stations  i  _>  1  and  j  =  1)  .   Assuming  that 
Ay.  .  =  An.  .  and  Ax.  .  =  As .    (i.e.  to  first  order  the 
mesh  is  rectangular) ,  and  that  the  streamline  shape  is 
known  explicitly,  the  slope  in  finite  difference  form  is 
expressed  as: 


dx 


. ,     _  nl+i,r1/2Anl+i,i  -  "i,i  +  ^2  K,i 


(B13a) 


or: 


1  ~  j-l  .  i  ~  I   ~  J"1  ~  ! 

oAn.'   ,    ,  +    Z  An.*   ,    ,  -oAn.'   ,    .  -^An!   ,   -     Z  An!    .  +  -=An!    . 

2  1+1,1     i=1     1+1,1     2     i+l,:     2     i,l       i=1     irj     2     1,3 

j         = 

1,3  As!    . 

lfD  (B13b) 


where 


nu  =  ki  =  ki  + 1}  {B14a) 


i      j_i       i  ~    j-1  ~        i 

n!    .  =  ±An!   ,  +   Z  An!   - ,,  =^An!   ,  +    Z  An!   ,  ,,  +  (j  -T) ;     i  >  1,    (B14b) 
1,3       2     i,l     -  1     i,l+l     2     i,l     ^_1     1,1+1  2  -  ^ 


The   inverse   of   the  radius   of   curvature    (v)    is   defined   in 
Eq.     (B15)    as: 
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v  -  5  =  dV**2   3/2  <B15> 

(1  +  (dy/dxp)  6//L 


By  nondimensionalizing  and  substituting  the  slope  obtained 
in  Eqs.  (B13a)  or  (B13b)  into  Eq.  (B15) ,  the  following 
finite  difference  equation  is  obtained: 

1,3    i-l,]' 


An,  .        As! | 

Vi    i    =    R1  =  °  2  3/2  {B16) 

x'3      R      (1  +  <<5ifj)V/2 

Because  the  magnitude  of  the  slope  6  .  .is  never  greater 

—  6 
than  10   ,  it  can  be  neglected  in  the  denominator  simplifying 

the  expression  for  the  curvature.   Substituting  Eq.  (B13b) 

into  Eq.  (B16)  and  retaining  only  terms  of  first  order,  the 

curvature  becomes 


v! 


i,D 


^AnV,   .+   Z  An.'   ,    ,  -iAn'        .-An?   ,-2     Z     An.'   , 
2     1+1,1     <^_1     i+l,l     2     1+1,3         i,l         ^=1       1,1 


!  .  J-l     .  i  ~ 

+  An!    .  +^An!   ,    ,  +    Z     An.'       ,  -^An!   .    . 
1,3     2     1-1,1     £=1       1-1,1     2     1-1,3 


(B17) 


/(As!    .)2 

i,3 


If  the  mesh  is  not  specified,  Eqs.  (B9)  through  (B12) 
represent  five  equations  in  five  unknowns  and  are  nonlinear 
since  An',  p'  and  u'  are  dependent  functions.   When  the 
streamtube  shapes  are  specified,  An1  and  the  curvature  v' 
are  known  quantities  and  the  equations  become  linear  in  p' 
and  u' .   Therefore,  Eqs.  (BlOa)  and  (BlOb)  represent  two 
algebraic  finite  difference  equations  solvable  for  p'  and  u' 
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throughout  the  flow  field,  provided  suitable  boundary  condi- 
tions are  known.   For  example,  if  the  values  of  p'  and  u' 
are  known  at  points  a  (i,j),  b  (i+l,j)  and  d  (i,j+l)  in  Fig. 
3,  the  pressure  and  velocity  perturbations  can  be  calculated 
algebraically  at  point  e  (i+l,j+l).   This  procedure  is  con- 
tinued ("marched")  throughout  the  mesh  until  all  flow  pro- 
perties have  been  determined. 

The  remaining  flow  properties,  density,  and  temperature 
are  then  calculated  from  Eqs .  (B9)  and  (B12) ,  the  continuity 
equation  and  equation  of  state.   Once  this  has  been  accom- 
plished, the  energy  equation,  Eq.  (Bll) ,  is  checked  to  see 
if  the  properties  calculated  from  the  specified  streamline 
shapes  agree  with  the  desired  heat  distribution.   If  it  does 
not  coincide,  an  iteration  procedure  is  performed  in  which  the 
streamline  shapes  are  varied  until  agreement  is  obtained. 

Boundary  conditions  are  needed  along  at  least  one  stream- 
line (usually  a  wall)  and  one  normal  (usually  freestream 
conditions) .   In  the  transonic  thermal  blooming  problem,  the 
boundary  values  are  chosen  as  the  unperturbed  freestream 
properties  upstream  of  the  heat  release  region  and  the 
velocity  or  pressure  on  a  bounding  streamtube  which  is  far 
enough  away  from  the  heat  release  region  that  it  can  be 
considered  isentropic,  Fig.  2.   With  the  pressure  and/or 
velocity  known  on  two  of  the  boundaries,  the  solution  for 
all  the  flow  properties  can  be  calculated  as  mentioned  pre- 
viously.  If  a  coundary  condition  is  known  at  other  points 
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in  the  flow  field,  the  resulting  properties  from  the 
"marching"  procedure  must  agree  with  it  or  the  specified 
streamtube  configuration  must  be  altered  until  coincidence 
between  the  boundary  conditions  and  heat  distribution  is 
obtained  throughout  the  flow  field. 
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APPENDIX  C 

DERIVATION  FOR  THE  RELATION  BETWEEN  AREA 
VARIATION  AND  HEAT  ADDITION  AT  PRECISELY 
MACH  1.0  FLOW 


1.   AREA  CHANGE  WITH  HEAT  ADDITION  AT  THE  SONIC  POINT 

From  the  influence  coefficients  for  constant  specific 
heat  and  molecular  weight  [Ref.  27],  the  change  in  Mach 
number  with  heat  addition  and  area  variation  is  given  by: 


4  -  f4 

l-M 


where 


G(X)  =  M2(l+i(Y-DM2)  [(1+yM2)   dx  °   -  2^lnA)  ]  ,   (Clb) 


d(lnT  )/dx  and  d(lnA)/dx  are  assumed  to  be  functions  of  x, 
the  streamwise  coordinate,  and  dx  is  always  positive.   Based 
on  the  relation  between  streamline  shape  and  total  temperature 
change  obtained  from  the  linear  subsonic  flow  and  supersonic 
flow,  Figs.  4  through  7,  the  following  relationship  between 
area  variation  and  total  temperature  change  (heat  addition) 
is  chosen: 


, ,,   .    ,  d(lnT  ) 

dx     =  ~(Y+D  (l+6(x))  ^ 2-  .  (C2) 
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3  is  a  variable  providing  an  additional  degree  of  freedom 
in  determining  the  relationship  at  Mach  1.0.   Substituting 
Eq.  (C2)  into  Eq.  (Clb) ,  G(x)  simplifies  to: 

o        1               oo                                d(ln  T  ) 
G(x)  =  NT(l+i(Y-l)*0  ((M-1)y  -  (Y+D3)  ^ — .     (C3) 

Since  the  denominator  of  Eq.  (Cla)  tends  to  zero  in 
the  limit  as  the  Mach  number  approaches  unity,  the  numerator 
must  go  to  zero  to  permit  a  continuous  passage  from  subsonic 
to  supersonic  speeds  (or  vice  versa) .   The  condition,  when 
G(x)  passes  through  zero  simultaneously  with  the  Mach  number 
becoming  unity,  is  designated  the  critical  point  G(x)  and, 
from  Eq.  (C3) ,  is  equal  to: 

1      2  d(ln  V 
lira  G(x)  E  G*(x)  =  -j(Y+D  Bgj —  -  0     (C4) 

M-*l 

Hence,  for  heat  addition  -d(ln  T  ) /dx  >  0,  3  must  equal 
zero  or,  equivalently,  the  area  change  is  directly  propor- 
tional to  the  heat  addition  (i.e.  change  in  total  tempera- 
ture) for  steady  flow  from  Mach  1.0,  Eq.  (C5)  is: 


a  /i   tv\    i        d(ln  T  ) 

|iiiLAL  .  i(Y+1)  ^^-.  ,c5) 


2 

At  the  critical  point,  dM  /dx  is  of  the  indeterminate 

2 

form  0/0.   Therefore,  the  limiting  value  of  dM  /dx  at  the 

critical  point  is  found  by  applying  L'Hospital's  Rule  to  the 
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right-hand   side  of  Eq.     (Cla) . 


lim^E    (<*£,*    =    (dG/dx)* 
M^1    ax  dx  ~(dM2/dx)* 


(C6a) 


orf  equivalently, 


<%*•-*/- 


(— )* 


(C6b) 


Eq.  (C6b)  suggests  that  (dM  /dx) *  can  be  positive  or  negative 
and  the  flow  can  become  either  subsonic  or  supersonic  from 
the  sonic  point.   The  requirements  on  $  follow  from  Eq.  (C6b) 
where  G(x)  is  given  by  Eq.  (C3) . 


(— )* 


,  .   dG 

lim  3— 

M-l  dx 

B=0 


i(Y+D 


-(Y+l)(§|)*  +  Y(g-)* 


d(ln  T  ) 
o 

dx 

(C7) 


Thus,  at  the  critical  point,  the  flow  may  be  continuous 
through  the  sonic  speed  provided  (dG/dx) *  is  negative  and 
is  not  zero.   For  heat  addition,  Eq.  (C7)  can  be  solved  for 
the  critical  value  of  (d3/dx)*,  Eq.  (C8) : 


^dx;     y  +  1 


vdx  ; 


(C8) 


and  is  graphed  in  Fig.  1C. 
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Allowable  Solutions  for 
Heat  Addition  d  In  (TQ)  /dx  >  0 


W 


'Slope  y/(Y+D 


Subsonic ' 


^dx  ' 


-^-  Supersonic 


Allowable  Solutions  for 
Cooling  d  In  (TQ)  /dx  <  0 


FIGURE  1C. 


3  Variation  versus  Mach  Number  Variation 
for  Sonic  Flow 


From  the  sonic  flow  condition,  g  must  be  positive  for 
supersonic  flow  and  can  be  either  positive  or  negative 
for  subsonic  flow.   The  change  in  3  with  Mach  number  squared 
must  be  greater  than  y/(y+1)  for  either  flow  condition, 
supersonic  or  subsonic. 

In  conclusion,  for  sonic  flow  the  area  variation  is 
directly  proportional  to  the  change  in  total  temperature, 
Eq.  (C5) ,  where  the  change  of  3  from  the  sonic  point  must 
satisfy  Eq.  (C8)  . 
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2.   STREAMLINE  CONFIGURATION  FOR  MACH  1.0  FLOW  WITH  A 
GAUSSIAN  HEAT  INPUT  DISTRIBUTION 

In  the  first  part  of  this  Appendix  the  streamtube  area 

change  was  shown  to  be  directly  proportional  to  the  change 

in  total  temperature  at  precisely  Ma  =  1.0;  see  Eq.  (C5) . 

To  derive  an  expression  for  the  change  in  total  temperature, 

a  slightly  different  form  of  the  energy  equation  is  used, 

Eq.  (C9), 


pu  |^  (CpTQ)  =  Q  =  la  (C9) 


in  terms  of  the  total  temperature  T  .   Assuming  a  constant 
Cp  and  a  Gaussian  heat  intensity  distribution  in  the  form 
given  in  Eq.  (2) ,  Eq.  (C9)  becomes: 

ni,j 

dTo.  -la  ,2  s.'2 

iti  =  2Q  g,       /   exp  -  -Z-y  dy'  exp  -  -^  ds! 


To,  .    fc~«"l,j     ,  20'" 

1,3  i,j 


(CIO) 


where  n!  .  represents  the  total  displacement  of  the  jth 
streamtube  from  the  centerline  (j=l) ,  and  QTO,  defined  by 
Eq.  (12),  is  evaluated  at  M^  =  1.0.. 

Equation  (CIO)  is  now  substituted  into  Eq.  (C9)  to  give 
an  expression  for  a  change  in  total  streamtube  growth  over 
a  streamline  distance  As?  at  M  =1.0. 

1         oo 

n!  .  2 

3n!  .    ,  ..._        x'3         .2  s! 

___Ll1  =  \l+1]la  (  exp  -  -*-_  dy  exp  -  -i,  ds! 

ni,j    4Q=°nl,j      J  2a'2  2a'2   X 

-n!  . 

1'3  (Cll) 
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or  more  conveniently  for  computational  purposes: 


An' 


dn!  . 
if  1 


i,h 


j     exp 


-_Y. 


,2 


20 


,2 


dy' 


"An!  . 


,2 


exp  - 


20 


-*■   ds! 
,2    i 


1        In°t  S- 

t(Y+1)  o2-  exp ±~  ds'   ;   i>l,  j=l 

H  ^oo        20 %  L 


(C12a) 


and 


dni,j  =  dni,i  +  fr+U   Jf! 


r  j-1    (2£+l)An.  £/2 

I        / 

1   -(2£-l) (An±  £/2) 


exp 


•  2 

-  ^—  dy 


20 


,2 


,2 


exp  - 


2a 


-*■  ds! 
,2    l 


i       I  a   j-1 
dn!  x  +  i(Y+l)  g2_   J   exp 

£=1 


(*An!  -) 


,2 


2a 


,2 


exp  - 


2a 
(C12b) 


r?dsi 


for  streamtubes  off  of  the  centerline  (j  >  1) .   Equations 
(C12a)  and  (C12b)  can  be  integrated  to  give  an  expression 
for  the  total  streamtube  perturbation  from  its  initial 
width  n'  .  as: 


An!  Ioay/2"  """   (       s! 


;  i  >  1 
(C13a) 
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i  j-1 

n!  .  -n'    =  T  An'    +   T    An! 


A£  I 


lAni/j  +  l(Y+l)  ^_ 


1=1  £2  sj 

L        exp 

£=1        2 


— -  (i  +  erf  — 1 — )   ;   i>l,  j>l  . 
o'Z  /2  a' 

(C13b) 


An  alternate  form  of  Eqs.  (C13a)  and  '(C13b)  which  can  be 
used  in  Eqs.  (B13a)  and  (B13b)  to  simplify  the  expressions 
for  the  slope  and  the  radius  of  curvature  of  the  streamlines 
is  as  follows: 

(Y+l)/2  a'  la  2  t  s!   ^ 

Afli.j ^ exP("  -^75-'  t1  +  erf  5^  }  ' 

i>l,  j>l        (C14) 

where  it  has  been  assumed  throughout  that  the  initial 

streamtube  widths  are  uniform;  that  is,  An,1  .  =  1  or 

J-  /  J 

An'  .  =  0  for  all  j . 

An  expression  for  a  streamtube,  which  is  far  enough  from 
the  heat  release  region  that  it  can  be  considered  isentropic, 
is  now  derived.   Carrying  out  the  integration  over  y'  in 
Eq.  (C12a)  and  assuming  that  n!  .  >  4a'  an  expression  for 
the  slope  of  the  bounding  streamtube  is  obtained: 


dn!  .*    (y+1)I  a/2lf  a1         s!  2 

-affJ- j2 exp  -  -±2-  ;  i>l,  j>l.    (CIS) 

i,  j*  °°  2a1 
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Integrating  Eq.  (C15)  with  respect  to  ds!  gives  an  explicit 
equation  for  the  bounding  streamtube  perturbation  from  its 

initial  position  n,'  .  *  as: 

1/] 


1  j*r_1 


(Y+DI   aa'2  s! 

j2 (i  +  erf-i-)       .  (C16) 

4y«>  {ia  » 
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APPENDIX  D 

INTEGRAL  EQUATION  DERIVATION  FOR  TREATING 
SMALL  DISTURBANCE  TRANSONIC  FLOW  WITH  SHOCKS 
FOR  AN  INFINITE  WALL  WITH  GAUSSIAN  SLOPE 


1.   GENERAL  DERIVATION  OF  INTEGRAL  EQUATION 

In  this  appendix,  an  integral  equation  is  derived  for 
the  small  disturbance  transonic  flow  with  shocks  on  an 
infinite  wall  with  Gaussian  slope.   The  small  perturbation, 
steady,  two-dimensional  transonic  equation  without  heat 
addition  can  be  written  as: 


(1"m2)*x'x'  +Vy«  =  (1-Mco2-M«2^+D<J)x,)c})x,x,  +<i>y.y.  =  0 

(Dl) 


where  the  local  Mach  number  is  approximated  by  the  freestream 

2 

Mach  number  plus  a  perturbation  quantity  (M   (y+1)  <$>    , )  . 

Equation  (Dl)  is  valid  only  in  regions  where  the  necessary 
derivatives  exist  and  are  continuous.   Therefore,  the  integral 
equation  derivation  includes  an  equation  for  the  transition 
through  a  shock;  the  equation  is  the  classical  relation  for 
the  shock  polar  [Ref.  34],  Eq.  (D2) . 

2 

u  uJ  -  a ' 

v« 2  =  (u:  -  u; ) 2  — 5-2 (D2) 

b      a   b    ,,,|2 

/us  2 

—4t   -  u'u'  +  a,z 
y  +  1    a  b 
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where  subscripts  a  and  b  refer  to  conditions  ahead  of  and 
behind  the  shock.   If  a  small  perturbation  analysis  is 
carried  out  on  Eq.  (D2)  similar  to  that  performed  in 
obtaining  Eq.  (Dl) ,  the  following  relation  is  found  for 
the  perturbation  velocity  components  across  the  shock  wave: 

(1-M2)  (u£  -  u£)2  +  (v^  -  v£)  2  =  M^2  (Y+D  (u^  +  u£)  (u^  -  u£)  2/2 

(D3) 

Equation  (D3)  corresponds  to  the  shock  polar  curve  for 
shock  waves  of  small  strength  that  are  inclined  at  any  angle 
between  that  of  normal  shock  waves  and  that  of  the  Mach  lines 

The  applicable  boundary  conditions  for  the  integral 
equation  analysis  are  that  the  perturbation  velocity  must 
vanish  far  upstream  of  the  disturbance,  Eq.  (D4) : 


3<fr 


94> 
'  3y' 


=   0  (D4) 

x=_oo 


and  that  the  flow  must  be  tangent  to  the  wall  surface, 
Eq.  (D5) : 

y.=0  1'3  Z° 


where  Sz'/^x*  is  the  slope  of  the  wall  surface  in  the  x'  (s') 
direction  which  can  be  approximated  by  Eq.  (C15) ;  and  x  is  an 
equivalent  thickness  ratio  for  the  surface  defined  as  the 
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maximum  wall  growth  t  (cm)  in  the  y-direction  divided  by  a 
characteristic  length  /2~¥  a  for  M  =  1.0. 

00 

2 

.  (-y-I)I    cnra 

t   =  - —      ;  t  =  - (D6) 

/2tt   a  'Poouoo 

Also,  it  is  necessary  to  prescribe  that  the  direct  influence 
of  a  disturbance  in  the  supersonic  region  proceeds  only  in 
the  downstream  direction. 

Sprieter  and  Alksne  [Ref.  10]  derive  in  detail  the 
integral  equation  by  a  Green's  function  analysis  for  transonic 
flow  around  bodies  having  subsonic  freestream  velocities,  and 
a  local  shock  discontinuity.   The  nonlinear  term  in  the 
differential  equation  was  maintained.   This  derivation  can  by 
extended  to  the  streamtube  configuration  of  the  boundary 
condition  for  the  internal  region.   For  the  sake  of  brevity, 
only  the  important  equations  and  pertinent  assumptions  will 
be  given  in  the  following  analysis. 

Since  the  object  of  the  analysis  is  to  determine  the 
pressure  perturbation  and/or  velocity  perturbation  on  the 
bounding  streamtube,  an  equivalent  equation  for  the  velocity 
perturbation  can  be  readily  obtained  by  differentiating 
Eq.  (Dl)  with  respect  to  x1;  it  is: 

2     ~  2 
(l-M„2)^,x,  +  8J,  .  =  VlJ(y+l)    1— j  (Hi-)        (D7) 
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For  the  Green's  function  analysis,  it  is  advantageous  to 
normalize  Eq.  (D7)  with  the  following  linear  transformation; 

x  E  x'     ;    y  =  By' 


(D8) 


M2    (Y+l)u» 

• 

M^y+Du' 

2 

3^ 

~                  3 
3J 

where      3  =   /l  -  M 

In  this  way,  Eq.  (D7)  reduces  to  the  following: 


32u    32u  _  a2-     32   ,u2  /r>ox 

3xz    By"2  3xz   z 


where  u  can  further  be  defined  by  the  approximation  to  the 

2      2 

local  Mach  number,  Eq.  (Dl)  ,  namely  M  =  M^  (1  +  (y+l)u),  as 

M  2(y+1)u'    M2  -  M  2 

—       00     ■  oo 

U  =  *-  =  j  (D10) 

1  -  M        1  -  M  ' 


00  00 


From  Eq.  (D10) ,  it  is  clear  that  for  a  subsonic  freestream 
Mach  number:  u  <  1  when  the  local  flow  is  subsonic,  u  =  1 
when  it  is  sonic,  and  u  >  1  when  it  is  supersonic. 

The  Green's  function  analysis  is  a  Neuman  problem  for 
finding  a  function  u  in  a  bounded  region  R,  bounded  by  a 
finite  number  of  closed  smooth  contours  C,  such  that  Eq.  (D9) 
is  satisfied  in  R,  u  and  its  first  partial  derivative  are 
continuous  in  R  (shocks  are  excluded  by  judicious  selection 
of  the  bounding  contour)  and  the  normal  derivative  3u/3y 
is  prescribed  on  the  boundary  of  R,  Fig.  ID. 
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(5,-(y0+n)) 


FIGURE  ID.   Integration  Region  for  Green's  Function  Analysis 
on  Flows  Past  Bounding  Streamtubes  at 
Transonic  Speeds 


Green's  theorem  states  that: 


/   (uVG«n-GVu*n)  dC  = 


//  (uV2G-GV2u)  dR     (DID 
R 


where  the  directional  derivative  on  the  left  side  are 
taken  along  the  normal,  drawn  inward,  to  the  curve  C. 
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The  Green's  function  G  is  constructed  so  that  it  satisfies 
the  Laplace  equation  v2G  =  6  U~^ ,T]-Y+Y Q)    with  3G/3n  =  0  on  C, 


G(£,n,x,y-yJ  =  £-   ln(^i)  =  ±   ln(R) 


(D12) 


where 


R  = 


/(x-c)2  +  (y-yQ-n) 


/(x-T)2  +  (y+y0+n)2 


and  f  and  r\   are  the  variables  of  integration  while  x  and  y 
are  the  coordinates  of  a  point  P.   With  the  Green's  function 
of  Eq.  (D12) ,  Eq.  (Dll)  becomes: 


2    —2 


u<*,y)  -  I J    G  £-_  (^-)  dndl 


=   |(GB-u^) 


-«   3n 


an 


3u 


9G, 


aj  -  f     (g  ^  -  u  ^) 


>=0 


3n 


3n 


dn 


x=x 


3u   -  3G 


+  f    (G  ^£  -  u 

sb    3n    3n 


3u   -  3G, 


dn   +  /  (G  ^  -  u  -2^)  dC   (D13) 


x=x. 


an 


an 


where  the  respective  integration  regions  and  contours  are 
designated  in  Fig.  ID.   Integrating  the  remaining  surface 
integral  by  parts  twice  with  respect  to  f ,  applying  the 
Green's  function  boundary  condition  on  C-,  ("n=0  or  y=y" '  )  , 
allowing  the  radius  r  of  the  curve  C~  go  to  infinity,  and 
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decomposing  the  line  integral  over  the  shock  wave  into 
components  parallel  to  the  axes  of  the  coordinate  system 
yields  Eq.  (D14) . 


u(x,y)  -  T  u2(x,y)  = 


oo  oo       3^G    —  — 


0  -°° 


ar 


-'{ 


GiL  (u-|u2)  -  (u4u2)  iS 

a^T    2        2    91 


Gi-(u-fe2)  -  (u-fe2)  i£ 
aT   z  2   al 


c  3u  .  —  3G 


an 


an  J  a 


n    au    -  8G 
G u  — 


_  an 


anj 


cos (n, n) 
cos(n,T) 


dn 


(D14) 


where  u_   represents  the  linear  wall  velocity  u  when  the 
wall  shape  is  specified. 


au 


dG 


77   -    r    n     du     —  dU 

-°°    an    an 


dX 


(D15) 


n=o 


Sprieter  and  Alksne  [Ref.10]  in  their  derivation  delineate 
the  advantages  for  integrating  the  surface  integral  over  R 
by  parts  as  was  done  in  arriving  at  Eq.  (D14) .   The  advantages 
are  important  enough  that  they  are  reiterated  here.   First, 
the  double  integral  of  Eq.  (D13)  shows  a  very  strong  influence 
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of  the  velocities  in  the  region  immediately  surrounding  the 
point  P  since  they  are  multiplied  by  ln(R).   This  influence 
is  largely  nullified  in  the  double  integral  of  Eq.  (D14) 
because  part  of  the  region  has  a  negative  influence  and  part 

has  a  positive  influence.   The  predominant  influence  in  the 

1  —2 
latter  case  is  furnished  by  the  term  y  u  standing  outside 

the  integral.   Next,  the  contribution  of  distant  regions 

is  diminished  in  importance  in  the  double  integral  of 

Eq.  (D14)  since  their  influence  varies  inversely  with  the 

square  of  the  distance,  rather  than  with  the  logarithm  in 

Eq.  (Dll) .   A  further  advantage  is  that  the  value  of  the 

double  integral  of  Eq.  (D14)  is  continuous  through  a  shock 

wave  rather  than  discontinuous  as  is  the  case  with  Eq.  (D13) , 

Lastly,  a  point  of  great  importance  in  the  approximate 

solution  arises  from  the  fact  that  the  integration  by  parts 

1  —2 

provides  extra  terms  (those  containing  y  u  )  in  the 

integrals  along  the  shock  surface  S  which  combine  with  those 
already  present  in  such  a  manner  that  the  contribution  of 
these  integrals  becomes  very  small  when  the  shock  wave 
approaches  a  normal  wave,  as  is  usually  the  case  at  high 
subsonic  speeds. 

Therefore,  u  must  satisfy  Eq.  (D14)  where  u^  is  given 
by  Eq.  (D15)  and  the  shock  relationship  Eq.  (D16) ,  which  is 
the  normalized  form  of  Eq.  (D3) . 

("a  "  "b>  2  +  (7a  "  7b>  2  "  k   (»a  +  %»  ("a  "  V  2     (D16» 
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In  the  case  of  a  specified  symmetric  wall  shape,  Eq.  (D14) 
can  be  considerably  simplified  as  it  permits  the  introduction 
of  the  boundary  conditions  which  are  specified  at  the  outset 
of  the  problem,  Eqs.  (D17a)  and  (D17b) ,  into  the  integral 
over  the  wall  surface. 


u(x,0)  =  uLw  =  KNOWN 


(D17a) 


and 


3u 

3n 


9v 
3T 


(3v} 

3T 


-ii 


-a2        -e 

t  —=7r   exp  -  — 2—. 


(D17b) 


where  Z  represents  the  reduced  ordinate  of  the  surface,  and 
T  the  reduced  thickness  ratio  given  by: 


T  = 


M„  (Y+1)T 

e3 


(D18) 


Substituting  the  relationships  of  Eqs.  (D17a)  and  (D17b)  into 
Eq.  (D15)  completely  determines  at  the  outset  of  the  analysis 
the  linear  wall  velocity  distribution: 


^Lw^O)  =_lim   (  ^/Jln(^)  JL§) 


y-*y, 


dl 


n=o 


_  1  rm   dZ   dl 
11   -»  dl   (x-I) 


(D19) 
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2.   SIMPLIFICATION  AND  APPROXIMATE  SOLUTION  OF  THE 
INTEGRAL  EQUATION  FOR  TRANSONIC  FLOW 


To  obtain  a  solution  to  Eq.  (D14) ,  some  further 
assumptions  must  be  made.   They  are:   (a)  all  shock  waves 
lie  in  a  plane  transverse  to  the  beam,  and  (b)  the  shock 
waves  are  normal  (i.e.  normal  to  the  local  flow  direction). 
These  two  assumptions  reduce  Eq.  (D16)  to  an  expression, 
Eq.  (D20) ,  which  permits  an  advantageous  introduction  into 
Eq.  (D14)  to  eliminate  the  integral  over  the  shock  S,  Eq.  (D21) 


,-  1  -2v     /-  1  -2X 
(u"2  u  >a  =  (U"2  u  >b 


(D20) 


2-    (u-i  u2)   =  1-  (u-i  u2), 

3T    2   a   H        b 


oooo  _2   2      r 
u(x,y)  =  \{x,y)    +   |  u2(x,y)  +  ^  /  /  ^-  2_  in (ji)  dndl 

— °°0     3  C       2 


(D21) 


Equation  (D21)  will  be  used  as  the  basis  for  obtaining 
the  velocity  distribution  on  the  bounding  streamtube. 
Approximate  solutions  of  this  equation  are  obtained  numerically 
by  starting  with  an  assumed  velocity  distribution  and 
iterating  until  convergence  is  obtained.   If  mixed  flow 
exists,  the  initial  velocity  distribution  must  include  a 
proper  discontinuity  complying  with  the  shock  relations  of 
Eq.  (D20) . 
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Due  to  the  double  integral  in  Eq.  (D21) ,  an  iteration 
process  would  be  cumbersome  unless  it  could  be  reduced  to 
a  single  integral  by  introducing  suitable  approximations. 
Oswatitsch  [Refs.  35  and  36]  has  shown  that  approximate 
knowledge  of  the  velocity  distribution  in  the  vicinity  of 
the  surface  can  be  expressed  in  terms  of  the  local  coor- 
dinate y,  the  ordinates  of  the  wall  shape  z"(x)  ,  and  the 

desired  but  unknown  velocity  distribution  u  (x,  0)  on  the 

w 

surface.   This  permits  the  reduction  of  the  double  integral 
to  a  single  integral. 

Oswatitsch  has  considered  such  an  approximate  relation 
in  which  the  velocity  u(x,y)  starts  from  the  value  u  (x,0) 
at  the  surface  with  an  initial  rate  of  change  given  by  the 

irrotationality  condition,  Eq.  (D22) ,  and  vanishing  at  large 

—2 

distances  as  1/y  . 


au 
9y 


w 


8v 
8x 


(D22) 
w 


The  resulting  relation  takes  the  following  form: 


u  (x,0) 
u(x,y)  =  - t  (D23) 

(i  +  (y-yQ)/b) 


where  b  is  a  function  of  x  satisfying  the  irrotationality 
condition  at  y  =  y  .   Solving  for  b  using  the  boundary 
conditions  of  Eqs .  (D17a)  and  (Dl7b)  yields: 
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2u  (x,0) 

b(x)  =  -  -J* «-  .  (D24) 

dzZ/dx^ 


Substituting  Eq.  (D23)  into  the  double  integral  of  Eq.  (D21) , 

integrating  with  respect  to  r\   and  setting  y  =  y   results  in 

a  simplified  integral  equation  for  calculating  u  (x,0). 

For  computational  purposes,  the  integral  equation  is  expressed 

as: 


uw(x,0)  =  uLw(x,0)  +  |u2(x,0)  -  jl  (D25) 


where 


»  u2(£,0)    j_-       _ 
1=    /  -2-g E(^-F^)  d£  (D26) 

—  00 


and  the  function  E  is: 


,2  .  „4. 


E(^i)  =  E(X)  =  — =-r[|-|X|  (5-10X^  +X^)  - 

b  (1+X2)5  2 

(1-10X2  +  5X4)ln|X|  -^-(1+X2)  (25-71X2-X4-X6)] 


Although  I  is  a  function  of  u   and  is  unknown,  it  is  informative 
to  rewrite  Eq.  (D25)  by  solving  for  u"w  in  terms  of  I  and  u^, 
thus: 

u  =  1  ±Jl  -  (2uLW-l)   =  1  ±Jl  -  L   ,   (D27) 
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where 


L  =  2uT   -  1 

Lw 


a  known  function  since  uT   is  known  once  the  wall  shape  is 

LW 

specified.   It  should  be  noted  from  Eq.  (D27)  that  the 

discriminant  must  always  be  positive  in  order  to  obtain  real 

values  for  u"  ,  thus  I  >  L.   Furthermore,  the  choice  of  the 
w         — 

plus  or  minus  sign  determines  whether  the  local  velocities 
are  subsonic  or  supersonic.   A  change  in  sign  at  a  point 
where  the  radical  is  zero  corresponds  to  a  smooth  transition 
through  the  sonic  speed.   A  change  in  sign  at  a  point  where 
the  radical  is  not  zero  corresponds  to  a  discontinuous  jump 
in  velocity,  namely,  a  normal  shock  from  supersonic  to  sub- 
sonic flow  when  progressing  in  the  flow  direction.   A  dis- 
continuity in  the  reverse  direction  is  inadmissable  since 
it  corresponds  to  an  expansion  shock,  an  impossible  phenomenon 
which  violates  the  second  law  of  thermodynamics. 

3.   COMPUTATIONAL  TECHNIQUE 

In  this  section,  a  brief  description  of  the  computational 
technique  used  to  solve  Eq.  (D25)  will  be  presented.   Sprieter 
and  Alksne  [Ref.  10]  describes  in  great  detail  how  a  judicious 
selection  of  the  initial  velocity  distribution,  depending 
upon  the  surface  shape  and  flow  configuration  (i.e.  entirely 
subsonic  flow  or  mixed  flow  with  possible  shocks) ,  is  essen- 
tial in  providing  a  rapid  convergence  to  the  solution.   It 
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will  be  assumed  for  the  following  discussion,  therefore, 
that  a  judiciously  selected  velocity  distribution  is  known. 
The  method  of  solution  is  as  follows: 

a.  From  the  known  surface  shape  in  normalized  coor- 
dinates Z,    calculate  its  first  two  derivatives  with  respect 

to  the  normalized  flow  coordinate  x;  that  is,  dZ/dx  and 

2—  —2 
d  Z/dx  .   Since  the  normalization  depends  upon  the  freestream 

Mach  number,  which  is  unknown,  the  value  of  7,  Eq.  (D18) , 

will  be  a  parameter  in  the  iteration  technique. 

b.  Calculate  the  linear  wall  velocity  u_  (x, 0)  from 

JjW 

Eq.  (D19) .   As  t  does  not  depend  on  Xr    it  can  be  taken 
outside  of  the  integral  of  Eq.  (D19)  giving  an  expression 
directly  proportional  to  the  parameter  7. 


uLw(x,0)  =  x  WLw(x,0)  (D28) 


where  W   (x, 0)  is  a  function  of  only  the  position  x  and  is 
fixed  once  the  surface  shape  is  specified.   The  quantity 
L  follows  immediately  from  Eq.  (D27) . 

c.  Assume  an  initial  value  of  t  consistent  with  the 
velocity  distribution  and  calculate  b(x)  from  Eq.  (D24) . 

d.  Since  b(x)  is  now  known  and  an  initial  wall  velocity 
u  (x,0)  has  been  selected,  Eq.  (D26)  can  be  integrated 
numerically  to  get  I. 

e.  The  I  function  calculated  for  the  initial  value  of 
7  is  graphed.   7  is  then  varied  in  the  L  function  until 
the  two  curves  are  tangent  at,  at  least,  one  point  (i.e. 
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I-L  =  0)  .   As  was  stated  earlier,  I  >  L  at  all  points  x 
along  the  surface  for  the  velocity  to  be  real. 

f.  Using  the  L  function  that  satisfies  the  tangency 
requirement,  a  new  approximation  to  the  solution  is  obtained 
from  Eq.  (D27) ,  where  the  sign  chosen  for  the  radical  is 
important. 

g.  The  iteration  procedure  continues  with  the  new 

u  (x, 0)  and  T  until  both  quantities  converge  to  a  solution. 

h.   Once  a  solution  is  obtained,  the  freestream  Mach 
number  for  a  given  flow  condition  is  calculated  from 
Eq.  (D18)  and  the  local  Mach  number  distribution  from 
Eq.  (D10) . 
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APPENDIX  E 
HODOGRAPH  TECHNIQUES  FOR  TRANSONIC  FLOW 

The  hodograph  plane  represents  a  coordinate  transfor- 
mation in  which  the  velocity  components  are  the  independent 
variables.   Spatial  coordinates  x  and  y  are  replaced  by 
velocity  components  u  and  v;  that  is,  the  physical  coordi- 
nates are  represented  by  x(u,v)  and  y(u,v)  [Ref.  26].   This 
coordinate  transformation  is  convenient  because  it  enables 
the  simple  presentation  of  data  or  solutions  and,  more 
importantly,  it  linearizes  the  two-dimensional  planar  tran- 
sonic flow  equation  that  is  nonlinear  in  the  physical  plane. 
This  equation  is  known  as  the  Tricomi  or  Tricomi-Euler  Equa- 
tion.  Boundary  conditions  are  difficult  to  satisfy  [Ref.  15] 

The  small  perturbation  transonic  equation  without  heat 
addition  can  be  written  as: 


<^>if^tf  -<(T+i)--ii;  w« 


where  for  irrotational  flow 


lil  -  221  .  o  (E2) 


Using  the  hodograph  transformation,  Eqs.  (El)  and  (E2) 
become: 
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(l-l£>2Zl  +  Ml  =  M2(Y+l)u»  -^1  (E3) 

3v'    3u'.  3v' 


and 


Ml  -  iZl  =  0  ,E4, 

3v»    3u' 


Eqs.  (E3)  and  (E4)  are  the  linear  transonic  hodograph 
equations. 

As  is  mentioned  in  the  introduction,  the  linearized 
small  perturbation  equations  for  supersonic  flow  are  hyper- 
bolic and  those  for  subsonic  flow  are  elliptic.   The 
transonic  hodograph  equations  change  from  elliptic  to 
hyperbolic  when  the  coefficient  of  3y'/3V  changes  sign: 


(1-M2)  -  M2(y+1)u'  =  0  (E5) 


The  critical  velocity,  where  the  change  of  sign  occurs, 
can  be  directly  solved  from  Eq.  (E5)  . 

(1-M2) 
u»*  = "  (E6) 

(Y+Dir 

Therefore,  for  u'  <  u'*,  Eq.  (E3)  is  elliptic,  and  when 
u'  >  u'*,  it  is  hyperbolic.   The  features  of  transonic  flow 
are  captured  by  the  transonic  hodograph  equations.   Since 
the  flow  occuring  in  the  external  region  and  on  the  bounding 
streamtube  is  considered  to  be  isentropic  two-dimensional 
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and  planar,  the  hodograph  method  gives  a  conceptual  idea  of 
what  the  resulting  flow  configuration  should  be  over  an 
infinite  wedge  of  prescribed  shape. 

Guderley  and  Yoshihara  [Ref.  15]  use  the  flow  angle  6 
and  n,  defined  in  Eq.  (E7)  ,  as  the  independent  variables. 


n  =  (Y+D1/3((w-w*)/w*)  (E7) 


/  2     2 
where  w  =  yu  +  v   .   With  this  selection  of  independent 

variables,  the  streamf unction  in  the  hodograph  plane  is: 


*nn  '  n  *ee  "  °  <E8» 


and  the  velocity  potential  function  is: 


*nn  "  "  *ee ■"  °  <E9> 


The  formulation  is  for  a  flow  with  constant  stagnation 
temperature  and,  hence,  can  be  applied  to  the  solution  for 
the  bounding  streamtube.   For  flow  over  a  wedge,  Guderley 
and  Yoshihara  found  a  steady  solution  for  M^  precisely  unity 

Fig.  lEa  shows  the  flow  over  the  bounding  streamtube 
in  the  physical  plane.   Points  A,  B,  C,  D,  and  E  are  shown 
mapped  in  the  hodograph  plane,  Fig.  lEb.   A  limiting  Mach 
line  is  shown  in  the  physical  plane.   Any  characteristic  in 
the  region  bounded  by  the  sonic  line  and  the  limiting  Mach 
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FIGURE  IE. 


Sonic  Flow  over  Bounding  Streamtubes 
(Physical  and  Kodograph  Plane 
Representations) 
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line  can  interact  with  the  sonic  line.   Any  characteristic 
downstream  of  the  limiting  Mach  line  does  not  interact  with 
the  sonic  line,  and,  hence,  the  flow  is  not  affected  by  the 
subsonic  flow  upstream  of  the  sonic  line. 
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APPENDIX  F 

DERIVATION  OF  THE  SLUING  RATES  NECESSARY  TO 
PRECLUDE  SIGNIFICANT  PHASE  DISTORTIONS  IN 
A  LASER  BEAM  AT  MACH  1.0 


When  the  phase  difference  between  two  points  in  a  laser 
beam  is  greater  than  a  quarter  of  a  wavelength,  the  deviation 
of  the  beam  from  its  original  direction  is  sufficient  to 
cause  considerable  distortion.   The  range  of  laser  sluing 
rates  necessary  to  prevent  such  phase  distortions  and,  hence, 
to  preclude  defocusing  of  the  beam  at  Mach  1.0  is  derived. 

Consider  a  laser  beam  that  is  slued  as  shown  in  Fig.  1. 
The  Mach  number  of  the  air  relative  to  the  beam  varies  along 
the  beam  as: 


H  =  \   =  2f   .  (Fl) 

a    a 


All  Mach  number  regimes  occur  along  the  beam.  If  L  repre- 
sents the  radial  distance  between  two  stations  on  the  beam 
where  the  Mach  numbers  are  M1  and  M2,  it  can  be  determined 
from  Eq.  (Fl)  that  L  is  equal  to: 


(M,-M,)a 

7f 


L  =    20  1         .  (F2) 


The  number  of  waves  (N)  of  laser  radiation  in  the  region 
of  length  L  is  given  by: 
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N  -  X  -  ^  -  <F3, 


where  A  is  the  wavelength  of  the  laser  beam  and  is  equal 
to  nv/c.   v  (radians  sec"  )  is  the  laser  frequency  equal 
to  approximately  2  x  10 13  Hertz  for  a  C02  laser,  and  c  (m 
sec   )  is  the  speed  of  light.   Therefore,  the  change  in  the 
number  of  waves  in  the  beam  length  L,  due  to  a  change  in 
the  index  of  refraction,  is  obtained  by  differentiating 
Eq.  (F3) . 


dN  =  vLdn 


(F4) 


The  change  in  the  refractive  index  is  related  to  the 
density  perturbations  in  the  beam,  a  known  quantity,  through 
the  Gladstone-Dale  Law  [Ref .  37] ,  a  simplified  version  of 
the  Lorentz-Lorenz  Law  for  gases  with  n  =  1,  Eq.  (F5) : 

n  =  1  +  kp  ,  (F5) 

where  k  is  the  Gladstone-Dale  constant  which  is  approximately 
equal  to  2.25  x  10   mykgm  for  air  with  wavelengths  beyond 
the  near  infrared  (X  >  ly) .   Differentiating  Eq.  (F5)  and 
substituting  it  and  Eq.  (F2)  into  Eq.  (F4) ,  the  desired 
expression  for  the  change  in  the  number  of  waves  is  obtained 
as : 


k(M0-M,  )avp    A 
DN  .  t_2_J: =_,&>  .  (F6) 

tic  p°° 
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It  is  of  interest  to  determine  the  range  of  sluing 
rates  for  which  the  change  in  the  number  of  waves  is  less 
than  a  quarter  of  a  wavelength  at  Mach  1.0  from  Eq.  (F7) : 


4k(M9-M.)avp 

°  ±  JT <^ax  -  5'in)  .        (F7) 


where  p *    is  the  maximum  positive  density  perturbation 
and  Pjjj^-  is  the  minimum  (maximum  negative)  density  perturba- 
tion in  the  beam  at  Mach  1.0,  points  A  and  B  in  Fig.  11, 
respectively.   As  an  example,  using  M~  =  1.001,  M..  =  0.999, 
the  flow  properties  in  Table  I,  and  the  maximum  difference 

in  the  density  perturbations  at  Mach  1.0  from  Fig.  11  of 

-4  -3 

1.5  x  10   ,  the  sluing  rate  must  be  greater  than  7.5  x  10 

radians/sec. 
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APPENDIX  G 

SUBSONIC  AND  SUPERSONIC  THERMAL  BLOOMING 
AND  COMPUTER  PROGRAM  (BLOOM) 


Tsien  and  Beilock  [Ref.  38]  have  obtained  solutions  for 
the  linearized  equations  of  motion  with  heat  addition  for 
subsonic  and  supersonic  flow.   The  solutions  are  for  a  line 
heat  source  q  in  a  uniform,  two-dimensional  planar  flow 
which  can  be  considered  to  be  a  perturbation  on  an  initially 
one-dimensional  flow. 

In  1971  and  1972  studies  were  conducted  at  NPS  by  Fuhs 
[Ref.  39]  on  the  applicatin  of  external  heat  to  the  reduction 
of  drag.   These  studies  required  thorough  analysis  of  heat 
addition  in  both  subsonic  and  supersonic  flows.   The  analysis 
of  the  flow  with  heat  addition  found  two  other  applications. 
One  application  was  external  burning  assisted  projectile.   This 
application  was  studied  by  CDR.  W.  J.  H.  Smithey  [Ref.  40]. 
The  other  application  was  recognized  by  Fuhs  [Ref.  4];  this 
involves  the  density  inhomogeneity  arising  from  quantum 
inefficiency.   Professors  Biblarz  and  Fuhs  [Ref.  5]  developed 
the  capability  and  extended  the  analysis  to  predict  density 
inhomogeneity  in  a  large  gas  dynamic  laser;  see  the  paper 
by  Fuhs,  Biblarz,  Cawthra  and  Campbell  [Ref.  6]  for  details 
of  the  experiment.   Subsequent  to  the  application  of  linear- 
ized solutions  to  flow  internal  to  a  laser,  it  became  apparent 
that  the  analysis  could  predict  flow  properties  in  thermal 
blooming.   In  a  series  of  papers  and  presentations,  Fuhs, 
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Biblarz,  Burden  and  Carey  [Refs.  7,  8  and  41-43]  developed 
the  appropriate  analysis  for  linearized  thermal  blooming. 

1.   SUBSONIC  FLOW  WITH  HEAT  ADDITION 

The  perturbations  to  the  flow  [Ref.  38]  are  given  by: 


a«  =       (Y-i)g  x 

2^V~^-    (X2  +  62y2}  <G1> 


v»    = 


(t-D  Sq 


r°°  °°      (x    +  3  y   ) 


(G2) 


P'         "    27ra    Bpm      J    2  ^  Q2    2.  (G3) 

00  ^°°      (x    +  6  y   ) 

(Y-l)MMq 

P'    = 3 2         2    2      "  3    *    6(y)     Kx)  (G4) 

where  the  symbols  are  defined  in  the  list  of  symbols  and 

/ 2~ 

3  =  /l  -  M   for  subsonic  Mach  numbers.   The  second  term  in 

Eq.  (G4)  is  the  wake  of  the  line  heat  source  q.   The 
temperature  perturbation  can  be  obtained  from  the  equation 
of  state. 

Biblarz  and  Fuhs  [Ref.  6]  in  applying  Eqs.  (Gl)  through 
(G4)  to  laser  internal  aerodynamics  have  presented  the 
integral  equations  used  to  solve  the  unbounded  heat  release 
region  problem  (i.e.  subsonic  thermal  blooming)  and  are 
listed  below  for  completeness.   Since  the  heat  source  is 
volume  distributed,  instead  of  a  single  infinitesimal  rod, 
an  integral  representation  was  used  to  account  for  all 
contributions  to  the  flow  properties  over  the  volume.   The 
geometry  involved  in  the  calculation  is  shown  in  Fig.  1G. 
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The  heat  quantity  h(x,y)  is  related  to  the  laser  output 
intensity  through  the  kinetic  lag  that  may  be  present.   To 
develop  an  equation  incorporating  the  lag  between  laser 
energy  and  heating,  it  is  necessary  to  examine  the  kinetics 
of  relaxation.   An  exponential  function  describes  the 
relaxation  with  a  relaxation  time  x [Refs.  44  and  45].   The 
heat  release  function  is  given  by: 


h(x,y)  =  la  J   f(x\y)  exp[-  (x"^'}]  dx'    (G9) 

—  00  00 


where  Q  =  la  f(x,y)   is  the  laser  energy  per  time  per 
volume.   The  volume  can  be  interpreted  as  an  area  times  the 
span  of  the  heat  addition  zone  across  the  planar  flow. 

If  significant  variations  of  beam  intensity  in  the  radial 
(beam)  direction  exist,  the  two-dimensional  heat  addition 
function  can  be  approximated  from  the  three-dimensional 
h(x,y,r)  by  an  average: 

ri+i 

h(x,y)  -  j£     I  h(x,y,r)  dr  (G10) 

ri 

where  Ar  =  r .   -  r .  .   The  beam  is  divided  into  segments  Ar 
in  length  and  Eqs.  (G5)  through  (G8)  are  solved  for  each 
segment. 

If  the  heat  release  h(x,y)  is  constant,  the  equations 
can  be  integrated  in  closed  form.   The  resulting  integrals 
are: 
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(Y-l)Ia 
G'    =   27TYP   3u     M<x'Y<*^n,B) 


(Gil) 


v'    = 


(Y-l)IQa 

2tiyp  gu    N<x'y<-^n,B) 


P'   =  - 


00^       00 


(Y-DM  I    a 

co     o 


2Tia  gp M(x,y;5,n,S) 

OO         •C^QO 


(G12) 


(G13) 


p'    =    - 


(Y-l)M  ±  a 

o 


3 M(x,y;£,n,B) 


(Y-Dla 

—5 //<5(y-n)   Kx-C)   dgdn 


aoo     UooPqo        £     f) 


(G14) 


where 

M(x,y;£,n,8)    = 


(y-n) 


In 


+    (-11  tan"l 


(x-g)2   +    82(y-n)2 
62(y-n)2 


g(y-n) 
(x-g) 


N(x,y;£,n,3)    = 


(x-g) 
2 
28 


In 


(x-C)2  +   82(y-n)2 

32 


n0  -i  C 


32  6 


(x-g) 
6 (y-n) 


The  first  term  in  Eqs.  (G8)  and  (G14)  correspond  to  the 
density  perturbations  in  regions  outside  the  wake  and  the 
last  term  represents  the  contribution  of  the  wake  and  can 
be  interpreted  as  the  consequence  of  the  heat  convected 
along  the  streamlines  by  the  flow. 
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2.   SUPERSONIC  FLOW  WITH  HEAT  ADDITION 

The  perturbations  to  the  flow  are  also  formulated  by 
Tsein  and  Beilock  [Ref.  38]  as: 


*'    =    "    2yp~6uq    5(X   -    By)  (G15) 

00  00 

*'    =   2yp^   Mx   "    W  <G16> 

(Y-D\q 

P'   "  2F-Hi 6(x  -   6y)  (G17) 

00  K  f  00 

(Y~1)M~q  (Y-l)a 

P'    = 3 6(x  -    3y)    -    U    XV* —  6(y)    I(x)  (G18) 


00  ^  00  00       oo       r  00 


Fuhs  [Refs.  4  and  39]  has  modified  these  equations  into 
an  integral  form  as: 


2 

a'   =  "  2yp"gu      J       h<5,n)    sin  v  ds  (G19> 

0 

s 

v«    =       (Y~1}       /        h(S,n)    sin   u   dS  (G20) 

2  ^PooUoo        o 


p'    =  ^ * /        h(S,n)    sin   y   dS  (G21) 

co^Poo 

0 
(Y-DMm        S 


p«    =  = /       h(£,n)    sin  y  dS 

2aB   Pp.      0 

00  00 


(Y-D 


Y~    J  I  MS,n)    «(y-.n)    Kx-5)   dSdn    (G22) 


Mmaoo      Poo  ^      n 

oo     oo         oo 


where  the  symbology  is  the  same  as  that  used  in  the  subsonic 
equations  except  that  B  =  /M2  -  1   for  supersonic  Mach  numbers 
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3.   COMPUTER  PROGRAM 

The  FORTRAN  computer  program  BLOOM,  included  at  the  end 
of  this  appendix,  was  written  to  compute  and/or  plot  the  flow 
perturbations  (u*,  v',  p'  and  p')  in  the  subsonic  and  super- 
sonic flow  with  heat  addition.  Using  a  modified  version  of  the 
computer  program  DIDER  [Ref.  5]  as  a  subroutine,  the  subsonic 
perturbation  equations,  Eqs .  (G5)  through  (G9) ,  and  the 
supersonic  perturbation  equations,  Eqs.  (G19)  through  (G22) , 
are  numerically  integrated  to  obtain  the  linear  thermal 
blooming  solutions.   The  program  calculates  the  two- 
dimensional  field  at  any  cross  section  of  the  laser  beam 
where  the  energy  release  distribution  is  known.   Finite 
kinetics  are  incorporated  as  described  by  Eq.  (G9) . 

The  program  (BLOOM)  is  versatile  in  that  it  can  calculate 
the  flow  quantities  for  several  subsonic  and/or  supersonic 
freestream  Mach  numbers,  heat  release  distributions,  relaxation 
rates,  and  freestream  flow  properties  in  a  single  computer  run. 
Since  the  length  of  time  for  each  individual  case  depends 
primarily  on  the  size  of  the  flow  region  and  on  the  mesh  that 
is  superimposed  upon  it,  time  requirements  must  be  determined 
on  an  individual  basis.   For  the  results  presented  in  Figs.  4 
through  7  and  Figs.  12  and  13  from  the  program  BLOOM,  the 
two-dimensional  region  of  interest  was  60  cm.  square  centered 
at  the  origin  of  the  Gaussian  heat  release  distribution  and 
was  descretized  into  a  mesh  of  cells  1  cm.  square  for 
subsonic  flow  (M  =  0.999)  and  into  rectangular  cells  whose 

00 

size  is  determined  by  the  Mach  angle  for  supersonic  flow 
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(M^  =  1.001).   Supersonically,  the  ratio  of  cell  height  Ay 
to  cell  length  Ax  is  equal  to  tan  u.   The  subsonic  results 
take  approximately  25  minutes  and  140  K  of  memory  (without 
plotting  subroutine  CONTUR)  while  the  supersonic  results 
take  approximately  3  minutes  and  285  K  of  memory  (without 
plotting  routine  CONTUR)  to  compile  and  execute  using  an 
IBM  360  computer. 

The  program  can  also  calculate  constant  heat  release 
distributions  for  single  and  multi-pass  laser  beam  configura- 
tions with  or  without  wall  reflections  and  kinetics.   In  its 
present  form,  the  program  can  not  calculate  supercritical 
(mixed)  flows  across  the  beam. 

A  main  feature  of  the  supersonic  portion  of  the  program 
is  that  it  sets  up  the  characteristic  lines  belonging  to  each 
individual  cell  (i.e.  their  diagonals)  and  adds  up  the 
contributions  for  each  flow  property  from  the  characteristics 
and,  in  the  case  of  the  density  perturbation,  from  the  wake 
at  each  cell.   The  calculations  along  the  characteristics 
are  facilitated  by  the  choice  of  cell  shape  and  size  as 
described  above. 
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